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ffcSeSDING PAGE BLANK NOT FILMED 

SECTION I 
INTRODUCTION 

This is an interim report concluding the third period of research 
under Grant NSG 1128, dealing with optimal design of subsystems of 

N 

the proposed Microwave Landing System. The research reported includes 
both optimal design studies of MLS angle-receivers jand a theoretical 
design-study of MLS DME-receivers . 

The angle-receiver results include an integration of the scan 
data processor and tracking filter components of the optimal receiver 
into a unified structure and then an extensive simulation study comparing 
the performance of the optimal and threshold (Phase III) receivers 
in a wide variety of representative dynamical and interference 
environments. The optimal receiver was generally superior, offering 
improvements ranging up to 20:1 or better in certain situations. 

The DME portion of this report includes a simulation study of 
the performance of the threshold and delay-and- comp are receivers in 
various signal environments. This study provides an -analysis of combined 
errors due to lateral reflections from vertical structures with small 
differential path delays, specular ground reflections with neglible 
differential path delays, and thermal noise in the receivers. 

The angle-receiver research and DME-receiver research are two, 
completely independent studies and are documented accordingly in the 
following eleven sections, the first part. Section's II thru VI, dealing 
xtfith the angle-receiver and the second part, Sections VII thru XII, 
dealing with the DME-receiver. No cross referencing occurs between 
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these major parts of this report, and also the symbol and notation 
systems used in the two parts are independent. For easy' access, 
however, results and conclusions from both parts of the year's work 
are summarized in Section XII. 
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PART ONE 


ANGLE-RECEIVER STUDY’ 
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SECTION II 

SIGNAL MODEL AMD APPROACH 

The reader is referred to our prior reports [1]/ [2], [3] 
for details . This is a summary included to communicate minor revisions 
in prior results and establish notation. 

GEOMETRY AND SIGNAL MODELING 

The angular coordinate to be estimated and other relevant 
quantities that evolve are assembled into an N-dimensional state 
vector x, modeled as the solution of a suitable- linear difference 
equation evolving in discrete- time, from scan-to-scan, and excited 
by a white zero-mean random process, {z(k), k = 0,1',*"'}. The receiver 
log-envelope signal, a continuous- time signal within a scan, is 
sampled throughout a window on each semi-scan at a sampling rate 
approximately equal to the i-f bandwidth and then suitably exponentiated 
and squared; the resulting J samples of the amplitude-squared 
envelope taken within a given scan are then normalized to a suitable 
measure of receiver -noise power and assembled into an observations, 
or measurement, vector u, which clearly is nonlinear in the state 
and also corrupted non-additively by receiver noise. For the kth scan, 
k = 0,1,2, , therefore we have the model form 

x(k+l) = F(k)x(k) + G(k)z(k) (2.1) 

u(k) = h(x(k), n(k)) 

relating state x, excitation z and observations u, generally. 

More specifically, in terms of a discrete-time variable local 
to the scan, and assuming the presence of a direct-path component. 
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a single multipath component and receiver noise, the jth component 

of u, say Uj , j = 1, , J, is approximated as follows with little 

error (see [ 


u. 

3 


{ap[6 - 0 A < T j)l + ct R pt® R “ ® A C T j)] cos (3 + “ sc Tj) + 

+ {a R p [0 R “ 0 A CXj ) 3 sin (p + w sc x..) + n g > 2 (2.2) 


where 

a - a (k) = -direct path signal-to-noise ratio 

0=6 (k) = angular coordinate of own A/ C 

a-r, ~ cio (k) = multipath signal-to-noise ratio 
K K 

0_ = 0_ (k) = angular coordinate of reflector specular point 
K K 


(2.3) 

(2.4) 

(2.5) 

( 2 . 6 ) 


6 = 3(k) = direct path-to-multipath phase difference at the 
beginning of the scan (propagating scan-to-scan 


as follows: 

3(k + 1) = 0 (k) + to sc T k 

where = time interval, kth-to-(k+l)th scans. (2.7) 
u = the scalloping rate (2.8) 

0^(. ) = the transmitting antenna scanning function (2.9) 

p [• ] = the transmitting antenna selectivity function (2.10) 

and 

n , n are independent Gaussian random variables 
c j s j with mean zero, variance 0.5. (2.11) 


On the basis of the above the state vector x is defined as follows: 
x = <a R ’ 0 R» 0 R» e »“s C ) T (2.12) 
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X % cl 

where ( ) denotes transpose and ( ) denotes ~ t - 

, dt 

Matrix F in (1) is then defined by 



where < > denotes mathematical expectation. This completes the 

modeling summary. 

APPROACH 

The objective of the desired MLS angle receiver is to produce 
an estimate of the A/C angular coordinate, 0, which is minimally 
affected by multipath interference. State estimation in conjunction with 
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the model (2.1) assumed is the approach used- to .develope the desired 
signal processor. Specifically, defining 


U' (k) = {u(k^), k^ = 0,1, ,k}, the sequence of observations 

from some initial time through the present; (2.16) 

A 

x (kjk) = estimate of x(k), given U(k) (2.17) 

Then the estimate evolution is described as follows : 

x(kjk) = x(k|k-l) + £(k|k) (2.18) 

where 


x(k|k-l) = F(k-l) i (k-l|k-l) (2.19) 

£(k|k) = estimate of the error in x(kjk-l), given U(k) (2.20) 

The calculation of £(k|k) in general, as defined, is complicated 
because of the form of h (’,*) in (1), as implied by (2). A 
simpler form for £ based in part on the "tracking assumption" that 
£ is "small" in some sense (and the vector Gw also) has been used 
thus far with good results (though without benefit of formal 
derivation from (2.20) as yet). This relation is 


£(k|k) = K(k) e (k | k) 


( 2 . 21 ) 


where 


e (k| k) = estimate of the error in x(k| k-1) in the neighborhood 

of zero error, given u(k). (2.22) 


A 

K(k) = a gain matrix, depending on x(k|k-l), Q.(k) , and 

statistics of £.(k| k) (2.23) 

The calculation of e (kj k) is based on the locally optimum 

estimation (LOE) criterion of Murphy [ 4 ] and does not make any use 

of the assumed dynamical model of state evolution, (2.1). This stage 
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of the computation processes the raw scan data {u(k), k = 0,1, } and 

is described in summary form in the next .section. The calculation 
of K(k) represents the determination of a suitable weighting matrix 
such that use of • the resulting £(k|k), (2.21), in the estimate update 
equation, (2.18) gives a smoothly evolving state estimate which .is 
appropriate to and consistent with the assumed dynamical model (2.1); 
in short, the constraints represented by the assumed state dynamics, 
heretofore ignored, are applied at this point. The quantity K(k) is 
the Kalman gain of a recursive tracking loop (or filter) , closed 
around the LOE. This aspect of the design is summarized in Section IV. 
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SECTION III 

SCAN DATA PROCESSOR DESIGN 

The observations taken on a 'scan have been modeled as a J-vector, 

u, with components u. , j = as given in (2.2). The application 

of locally optimum estimation to this problem requires the conceptualization 

of a noise-free observation vector, say q(k), whose jth- component 

q . has the form 
J 


A _ 

q . = u . 
'3 


f“pj(9) + « R Pj (6 R ) cos g^ ] 2 + [« R Pj (0 R )sin g..] 2 (3.1) 


n - 0 
c . 

3 

= 0 
n 

s . 

3 


= « 2 P j 2 (9) +2aa R p.(0) Pj (6 R ) cos + a R 2 P-j 2 ( 0 R ) 0.2) 


where 


g. ^ B + “sc T j 
J 


P j (0) = P [9 - 0 a (t )] 


(3.3) 

(3.4) 


and similarly for p^ (0 R ) . Using the q_. - formulization, it is 


possible to write u.. as follows: 


u . = q . +2 
J J 


q [n cos g. + n sing.] + n 2 4- n : 
J i J S 4 1 c s. 

J 3 3 3 


(3.5) 


= q. 


/ n. \ 

h! 

+ (n n ) 

c. S. j 

f n c. \ 

3 3 

2 

n 

\ S 3‘i 


l n s. 1 
3 


(3.6) 
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Letting the noise vector, (n n ) , be modeled as follows 

• c , s , 

3 3 


cos £. -sin 3. 
3 3 

sin cos 


where 


<S> = <S> =0 
{\ \) ' °’ v U1 
(vO - <v-i> ■* v 


clearly assures these same assumed properties in (n n ) , and 

c » s » 

3 3 

in addition simplifies u_. in (3.6) 


u . = q . + 2 n 
3 3 c . 

J 3 


+ n. 2 + n 2 
c . s . 

3 3 


Henceforth, references to receiver noise will refer to the new vector 
T 

(n H ) unless otherwise noted, 
c. » s . 

3 3 

Paralleling the prior development somewhat [ 3 ] , the likelihood 
ratio X(u|q), in terms of the revised formulations of u and q, is . 
given now by 


X(ujq) = H M (4q.u.) exp(-q ) 

j=l l J J J „ 


where, as before, M ( ) [for positive argument — as is the care here] 


12 . 



is defined in relation to I ( ) , the modified Bessel function of 

o 3 

the first kind, zeroth order, as follows 


M o (z) = I q (\/7“ ), z>0, real 


(3.13) 


The theory of locally optimum estimation is applied to the 
scan data processing problem by first assembling a selected subset of 
the state variables into a parameter vector Y and then processing 
the observations u to obtain an estimate, e, of the error in the 
current y-estimate, as follows: 

e = § ^ A (u|q) (3.14) 


where 


<P k 


<A (ujq) A T (u|i) 


(3.15) 



Y = Y 


(q(y), given in (3.2) 


(3.16) 


and 


A (ujq) = 


-f- In ’ 7(u | q(y) ) 

dy 


Y=Y 


(3.17) 


The estimate e is both locally unbiased at zero error (i.e. when 
Y = y) and has minimum mean square error of all estimates locally 
unbiased at zero error. The matrix § ^ is the \covariance matrix 
for the error in the estimate e when y = y. The averaging done in . 

rs 

the calculation of (3.15), is taken under the assumption that y- = y, 

A 

i.e. that the parameter value entering the q calculation equals the 
true value giving the u observation; the result is independent of the 

A 

observation and is a function only of q. The observations u, enter- 
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calculation of e only through A(u|q), which is a function of both 

A 

u and q. Substituting (3.12) into (3.17) we may write, specifically 


J 

I 

j=i 


J 

-I 


8q i 

9y 


3=1 

= D (q) w (u|q) 


where 


[ 


D(q) 


9q 1 3q, 


9y * 9y 


M 

w (u | q ) = 4u. (4q..Uj)-l 


[ o (4 V:l 

) -ij 

(3.18) 

\ 

4u. — — 

3 M 
J o 


(3.19) 



(3.20) 

t 

8q T ' 

97" . 

| , a matrix 

(3.21) 

.u.)-l 
3 3 

, a J-vector 

(3.22) 


and, as before 

h M ° (z) - 


(3.23) 


Using the previously developed asymptotic approximation for M^/M , 
i.e. 13, egn (III-42p)] 

M, 


1 / s ~ 1 

(z) .. 


M 




(3.24) 


The vector w(u|q), through which the : observations , u > enter the 
computations, has representative element 


w ■ 


(n|q) = 


u. 

3 


yl + q.u. 
3 3 


- 1 


(3.25) 
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Substituting (3.20) into (3.15) results in an expression for matrix 

A 

$(q), as follows: 

<Kq) = D(q) ^(u|q) w ' T (u|q) [ q^> D T (q) (3.26) 


D(q) H (q) D (q) 
w 


(3.27) 


where 


H w (q) = ^u|q) w T (u|q) |q^> 


(3.28) 


A simulation study was performed of the statistics of w (uj q) , (3 . 25) , 
involving up to 10,000 samples. Refer to Appendix A for details. 
Conclusions reached are as follows: 

- (u| q)|q^) = 0, independent of q (3.2 , 9 ) 

^(ujq) w fc (u|.q)fqy> - (q) (3.30) 


where 6., is the Kronecker delta, and 


h j <q) ‘ 


d^Culq) [ 


(3.31) 


The process (w^(u|q), j=l, J} is white; consequently. 


H^Cq) = Diag [^(q), h 2 (q), hj(q)J 


(3.32) 


The sample statistics study also produced an asymptotic approximation 

A 

for hj (q) : 


Y q > 


1 + 2 q . 

3 


(3.33) 
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whose error peaked briefly at 20% (for q ^ = 2) , but which produced 
good results in the receiver simulation studies. 

All receiver simulation studies to date dealt with a receiver 
design based on the following choice of the parameter vector Y : 


A 

Y = 


R 

5 R 


8 

The corresponding D-matrix is the following: 

2 


D(q) = 


(3.34) 


• ..2 a p j (0) + 2a R Pj (6>Pj (e R )cos^ . . . 

• • 

- • * 2 a 2 'P-j (0)pj (9) + (0)p.. (0 R )cos3.. . . . 

. . .2apj (e) P j (6 R ) cos 3^ + 2a R Pj 2 (0 R )... 

. . .2 aaR p^. (0 R ) pC 6 r ) cos + 2a R 2 Pj (0 R )p^ (0 R ) 

r ..-2aa p (0)p (0 ) sin 6 .... 

R J J R 1 


(3.35) 


where ( ) denotes d_ ( ). 

dt 
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SECTION IV 


TRACKING LOOP DESIGN 

A 

The scan data processor produces an estimate e of the error e 

A 

in the last . estimate y of parameter vector y, a masked version of the 

A 

state x. The algorithm employed provides that e is locally optimum 
at e = 0 in a least mean square error sense and supplies also the 
associated covariance matrix 

R = ^ (e r e) (e - e) T q ^ " y § 1 (q) (4.1) 

e=o 

The tracking loop design was obtained as follows: 

1. Generation of a pre-estimate of y by simply adding the 

estimate e to the value of y used in the calculation of e. 

2. Interpretation of the pre— estimate as a synthetic measurement 
of y , additively corrupted by white, discrete-time noise 
with covariance R. 

3. Use of a linear Kalman filter designed to accept the 
synthetic measurements produced (and the matrix R) , for 

A 

generating an update estimated, x of the full state. 

The method embodies certain assumptions, including specifically 

1. The whiteness of the synthetic measurement noise. 

2. A knowledge of matrix Q, defined in (2.15) 

3. The form of state-estimate error in (2.21). 

The algoritum therefore is deemed suboptimal; nevertheless the 
tracking performance results were good and probably represent the 
limiting performance obtainable from the standpoint ,of algorithmic 
complexity. 
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We will now formalize the approach mathematically and present 
in summary form the algoritum, which benefits from some simplifications . 


that are possible. 


Pre-estimate: 


where 


A * A A A ^ 

y(k|k) = y(kjk-l) + e(k) = Hx(k|k-1) + e(kj ' 

= [y(k)’~ e(k)3 + e(k) 

= Hx(k) + v (k) 

H = masking matrix associated with choice . 
of y 

v(k) = e(k) - e(k) 


and 


( T < k >) U 

<^(k)v T (k)^) | 


= 0 (locally unbiasedness) 

- R(k) = $ ”^(q (k) (by 4-1) 

e=o 


Kalman filter : x(k|k-l) = F(k-l) x (k-l[k-l) 

P (k|k-l) = F(k-l)P(k-l|k-l)F T (k-l)+G(k-l)Q(k-l)G' 
K(k) = P(k|k-l)H T [HP(k|k-l)H T + R(k)] _1 
x(kjk) = x(k|k-l) + K(k)[y(k[k) - H x(k|k-l)] 
P(k|k) = P(k|k-1) - K(k) HP(kJk-l) 

= [I-K(k)H]P(k|k-l) [X-K(k)H] T + K(k)R(k)K T (k) 


(4.14) is preferred to (4.13) for preservation of symmetry and 

\ 

positive-definiteness properties . 

Substantial simplification follows from ’substitution of (4.2) 
into (4.12 ) , giving 

x(k|k) = x(k|k-l) + K(k) e(k) 
where 


(4.2) 

(4.3) 

(4 .'4) 

(4.5) 

(4.6) 

(4.7) 

(4.8) 

(4.9) 

(k -1 ) (4.10) 
(4. -11) 

(4.12) 

(4.13) 

(4.14) 


(4.15) 
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(4.16) 


*K(k)e(k) = P (k| k-l)IT 5 ' [HP (k| k-l)H^ + $“ 1 T 1 <T 1 A = K® _1 A 
= PH T {$[HPH T + $" 1 ]}: 1 A 
- M(k)A (u | q) 

in whicli 

M(k) = P(kJk-l)K T [I + $ HPH 1 ]” 1 
Thus the state estimate update - operation, i.e. 

x(k|k) = x(k|k~l) + M(k)A(u|'q) , 
does not require the inversion of matrix $. 

The error covaviance update equation (4.14) can be written 
in terms of ® also by first noting from -(4.16) and (4.18) 
above that M = K® \ or more specifically 
K(k) = M(k) * 

^ -j 

Substituting this into (4.14) and recalling’ that R = ® ■ , gives 
P(k|k) - [I-MSH] P(k|k-1) [I-M®H] T + M3>M T 


(4.17) 

(4.18) 

(4.19) 

(4.20) 

(4.21) 

(4.22) 
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SECTION V 


SIMULATION STUDIES 

The principal simulation models with which we were concerned 

during the past year are the following: 

i. The environment and baseband receiver signal; 

ii. The LOE/Kalman filter recursive receiver structure, 

and specifically both multipath-adaptive. and non-adaptive 
variants, thereof; 

iii. A representation of the phase III MLS receiver, denoted 
the threshold receiver 

Simulation studies conducted, included principally the following: 

i. Crossing multipath interference, initiating as out-of- 
beam interference. 

ii. Time- varying in-beam multipath interference 

iii. Simulated landing scenarios 
Results are presented below and discussed; programs developed and 
used will be transmitted to the sponsor under separate cover. 
Simulation Models 

A. Environment and Baseband Receiver Signal 

Generally, the environmental dynamics are simulated with a 
state model of the form (2.1) (without the random excitation), using 
the state vector (2.12). To provide some commonality between the 
optimal and threshold receiver simulations, however, .the observations 
are generated in absolute-amplitude form. The full model is as 
follows : 

x(k+l) = Fx(k) , x(0) = x^ (5.1) 

v(k) = h v (x(k) ,a,n(k)) (5.2) 

where x = the initial state at the start of the simulation, 
o 
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F = the matrix (2.13) except T^ is a selected constant in 
the simulation (1/13.5 sec. for AZ, 1/40.5 sec. for EL). 

a = rms value of receiver noise at a point in the I-F 

channel having the same signal amplitude as the demodulator 
output. The parameter or is assumed known, being a receiver 
characteristic. 


V 


) = a matrix-valued function of its arguments which compiles 
the J-vector v(k) as one with representative element V j (k) , 
3=1, . , . , J, 
where . 


v. (k) 

3 



and u. 
J 


is as given in (2.2). 


All other quantities are as previously defined. The components of 
x q are sp'ecified either in program DATA statements or read in at 
run time . 


The quantity 6.., (3.3) is reduced to its principal value on 
(-31, II) after each change. 

Signal data samples are generated only during sampling windows 
of J/2 samples each, located in the TO and FRO scans respectively, 
and centered where the centroid of the received signal pulses are 
expected. For all runs to date 

J = 130 (5.3) 

corresponding (at the sampling rate of 160 kHz) to window widths of 8° 
in each semi-scan ^ 

B. The Optimal Receiver Simulation 

The optimal receiver simulation consists basically of the following: 
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i. Extrapolation of x to the present, via (4.9). 

ii. Scan data processor calculations of A, via (3.20), and 
$, via (3.27). 

iii. Kalman filter calculations as follows: 

(a) P(k|k-1), via (4.10) 

(b) Gain matrix, M(k) , via (4.19) 

(c) x(k|k), via (4.20) 

(d) P(k|k), via. (4.22) 

The scan data processor calculations begin with a computation of 
the squared amplitude observations vector u component-wise as follows 


v . 

u. --J- 
J 2a 2 

In the subsequent calculations the following models for the antenna 
selectivity function, p(0), and its derivative p(0) were used: 

• 0 = B/2.4 


r a 

4 


P(6) = / 


cos (1.2)II9/B 
^ l-(2.40/B) 2 


(5.5) 

elsewhere 


and 


( 


p(e) = \ 


+ 


V. 


0.3II 

B 


o.3n 2 


Signum (0), 

. r sin II (Z+l) 

It 2 

COS y (Z+l) - ■= 

f ( Z +D 


B 


cos ^ (Z— 1) 


1 <z+ « 

sin j (Z-l) 

' I (z-D 

2 


n 

2 


(z-l) 


2.49 

B 


0 = B/2.4 


(5.6) 


elsewhere 
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in which B, the -3db beam width in degrees, was given the value of 1 
degree. 

In the Kalman filter calculations the diagonal elements of the 
diagonal matrix Q were given values as follows: 

Qll = Q 3 g = max [0.25, 0.01a 2 (k) ,] (5.7) 


where a is the estimated direct path signal-to-noise ratio, thus giving 
some adaptation on the basis of 10% uncertainty in a (k + 1) , given a(k) 
(for a(k)>5) . 


Q99 - Q A a‘- |e | T # , 


(5.8) 


4 22 ^44 ' max' 

where T is the interscan. interval. All runs to date were made with the 

AZ receiver simulation, and, based on a prior study of representative 

landing patterns [3,p.40]., a value of 1 0 [ = 0.1 sec was used. 

' max 1 

JZ 


Q 55 - 0.04/T 


(5.9) 


corresponding to a mean-square uncertainty in oo sc (k + 1), given w gc Ck)> 
no greater than that which would cause an interscan extrapolation error 
in g of about 10 while tracking. 

The optimal receiver simulation is programmed with maximum 
dimensions of 8 and 5 for the vectors x and y respectively (and all 
associated matrices) . The actual dimensions used in the calculations 
however are parameterized with the integer variables NS and NG 
respectively. When NS=8 and NG=5, the optimal multipath-adaptive 
receiver, which has been described, obtains; when NS=3, NG=2, a lesser 
dimensional model of the same basic structure results for which 
x = (ct,0,0) T 
y = (a,0) T 
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(5.11) 



corresponding to an optimal receiver design predicated on the assumption 
of a multipath-free environment. Clearly, in the current work this 
latter design is a suboptimal one, nevertheless it served as a 
comparison standard in this work and is referred to as the suboptimal 
design or the nonadaptive design. 

C. The Threshold Receiver 


The threshold receiver simulation f.irst computes the log-amplitude 


envelope observations signal, v ^ og > component-wise as follows: 

V log = 2 0 l°g l0 (1 + Av ) 

3 


(5.12) 


where A = 100, corresponding to 40 db of logging action. The result is 
then filtered by a 25 kHz bandwidth low pass digital filter with transfer 
function 

H 25 (z) = 0-34831(1 + z" 1 ) (5 ,13) 

1 - 0.30336z~ 1 

The signal that results is then passed to the thresholding and 
interference-rejection logic that characterizes the standard phase III 
MLS angle receiver 'design. This is described as follows:' 

1. On each the TO- and FRO-semi-scans the signal peak within the 
tracking gate (located as described below) is determined 

and a threshold level 3db below the peak established. 

2. Dwell gates are defined for those intervals during which the 
log video signal exceeds the thresholds. The tracking gate for 
the next TO-scan will be symmetrically located about the expected 
dwell gate centroid position with a duration of 1.5 times the 
present TO-scan dwell gate duration; similarly for the FRO-scan. 

3. Dwell gates less than 15ys or greater than 350ys cause the scan 
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to be aborted and the prior estimate of 6 to be used again. 

4. When the dwell gates are acceptable, the interval from the center 
of the TO dwell gate to that of the FRO dwell gate is determined, 
quantized to 0.5ysec, and used to calculate the new estimate 0, 
which is output. 

The threshold receiver simulation was programmed to be interchangeable 
with the optimal receiver simulation as far as the main simulation program 
(and sampling window positioning) is concerned. Performance evaluation, 
however, was based on the angle estimate error filtered as follows: 

AZ: H(z) - — (5.14) 

1 - 0.5z 

EL: H(z) = ~' Z5 T~ (5.15) 

1-0. 75 z 

corresponding to available evaluation data. In this respect, however, 
the simulated threshold receiver is more like the phase II model than 
the phase III model (which apparently uses an a-0 filtered error for 
performance evaluation) . 

Simulation Runs and Results 

Four key parameters important to the performance of an MLS receiver 
are the following: 

S/N = Direct-path signal-to-noise ratio (denoted DSNRDB in the 

simulation), (db) . (5.16) 

p = Multipath-to-direct path signal amplitude ratio (5.17) 

F gc = scalloping frequency (Hertz) (5.18) 

0 = 0 - 0 p , the separation angle of the multipath interference (5.19) 

S6p iv 

The MLS receivers are expected to operate with S/N ratios, of 8 db or higher; 
values in the range 8 to 20 db were used in the simulation study. 
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A recent MIT study of multipath interference at air terminal areas [5] 
shows relative multipath amplitudes, p , to 1.0 or more and scalloping 


frequencies, F , to 1300 Hz. Representative values are as follows: 
s c 


p = 0.9, F = 2,22,51 Hz 
sc 

p = 0.8, F sc = 63,81,130 Hz 

p = 0.5, F = 381 Hz 
* sc 


(5.20) 


Values spanning these ranges of p, F gc were used in the study. Separation 
angles, 0 g , corresponding to both in-beam and out-of-beam multipath 
were considered. 

A fifth parameter, g, the r-f phase difference initial values 
(at the start of the simulation run) also affected results somewhat. Its 
effect is studied some and its value is always noted. 

A. Crossing Multipath Studies 

This scenario begins with 

0 - - 2.7 5° (AZ) 

sep 


de 


sep . n ,o, 

— r~ = + 0.7 /sec constant 
dt 


(5.21) 


and runs for 100 scans (approx. 7.4 sec). Forty runs each for the threshold, 
optimal and suboptimal receivers were made in this series, corresponding to 
various values of key parameters S/N, p, and F gc - In all runs the receivers 
were initialized in the track mode, i.e. all estimated variables produced by 
each receiver were initialized to true values. Figures 1 thru 10 show compara- 
tive results of selected runs for the optimal and threshold receivers; Fig. 1 
presents time histories of error for S/N = 14 db, no multipath; note the two 
plots are made to the same scale for easy comparison. Fig. 2 is the same, but 
with heavy multipath interference, p=0.8, = 51.3 Hz. (This scalloping 

rate and the associated value of 8, - 168°, produces maximum enhancement 
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Figure 1. Receiver Error Time Histories , Crossing Multipath 
Simulation, No Multipath (reference) case. 
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6 = 32.75° 


<? R = 31 .37° 


& R = 30.18° 


Figure 2. Receiver Error Time Histories grossing llultipath Simulation, 


S/N = 14 db, p 


.3, F gc = 51.3 Hz., 3 


-168.0 
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by the multipath on the TO scan and maximum cancellation on the FRO 

scan as the separation angle traverses zero). Figures 3, 4 and 5 

summarize for F = 5, 51.3 and 500 Hz, respectively, comparative 
s c 

studies of peak absolute error versus S/N and p. Figures 6, 7, 8 and 9 
present time histories of error for runs corresponding to selected' 
points in Figure 3, 4 and 5 representing both moderate and heavy 
multipath interference. The high ratios of improvement provided by the 
optimal receivers are especially noteworthy — typically about 20:1 
for the p=0.8 Cases. Figure 10 shows, for the S/N = 20 db, p=0.8 cases 
only optimal receiver results, plotted with enlarged scales to show 
more clearly the sample functions of the error process, which appears 
to be more random than that of the thresholds receiver (See 
Figures 6, 7 and 8). 

Tables 1, 2, 3, and 4 summarize all the crossing multipath studies 
Table 1 presents a comparison of the 3 receivers, optimal, threshold 
and sub-optimal (or non-adaptive) , in terms of peak-absolute-error. 
Tables 2, 3 and 4 each summarize the performance of respectively one 
receiver in terms of several error measures computed over the set of 
100 scans per run for each case. 

The crossing multipath scenario represents a strenuous test of the 
tracking capabilities of a receiver algorithm. Conclusions dram are 
as follows: 

1. The optimal receiver generally outperformed the threshold 
receiver, sometimes by a wide margin. 

2. The optimal receiver is algorethmicly much more complex than 
the threshold receiver. 
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Figure 10. Optimal Receiver Error Time Histories grossing 
Multipath Simulation ,S/N = 20 db, p = 0.3. 
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COMPARISON OF PEAK-ABSOLUTE-ERROR PERFORMANCE 
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{Table 2. Optimal Receiver Error Performance; 
Crossing Multipath Scenarios. 
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£able 3. Threshold Receiver Error Performance; 
Crossing Multipath Scenarios. 
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Sable 4. Non-adaptive Receiver Error Performance; 
Crossing Multipath Scenarios. 



4. Generally the suboptimal (non-adaptive) receiver performed as 
well or better than the threshold receiver. 

5. Generally the suboptimal design performes less well that the 
optimal design, though it was never observed to lose track. 

It is felt that the optimal and suboptimal receivers represent, 
within this family of receiver structures, two extremes in use of 
any information in the received signal concerning the multipath 
interference — both being generally superior to the threshold 
receiver. Further, it is felt that a carefully drawn design 
intermediate to these extremes can effect a substantial reduction in 
complexity at little loss in tracking performance with respect to the 
optimal design described. This, in part, is our recommendation for 
future work; the reader is referred to Section III for more details. 

E. Time-Varying In-Beam Multipath Scenario 

This is primarily a test of the multipath-acquisition capabilities 
of the optimal receiver. Figure 11 presents the result of a simulation 
run, which began with no multipath, the receivers tracking and 
S/N = 20 db. After about a second, multipath interference begins 
to appear at a constant separation angle of 0.5°, growing in amplitude 
to about 0.8 of the direct path signal amplitude, then diminishing, 
again to zero. The optimal receiver offers a 10:1 improvement in 
peak error performance. 
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Figure 11. Time- varying In-Beam Hultipath Scenario, 
Receiver Error Time Histories, 

S/N = 20 db, F = 5 Hz., g = 0°, 0 








C.’ Simulated Landing Scenarios 

-A landing scenario in general is characterised, in terms of our 
simulation' parameters, by simultaneously varying p, 0 g and ’ F gc [5] 

A case that was simulated is shown in Figure 12, suggestive ;of heavy 
in-beam multipath, a Fresnel reflection pattern and closing ranges. 
Error time-histories for the optimal and threshold receivers, operating, 
at 20 db S/N ratios, are shown in Figure 13. The receiver simulations 
were initialized in the track mode. 
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Figure 12. Representative Landing Scenario 


46 







SECTION VI 


THE EXPERIMENTAL SYSTEM 

A philosophy and design plan for an experimental system was 
described in the last report [3] . Three small circuit boards under 
construction at that time have been completed, but aside from that, 
no effort during the current grant period was expended on the experimental 
system. This task was halted because it was becoming increasingly 
apparent the computational requirements of the optimal receiver 
algorithm, as it was evolving, would he beyond the capacitates of 
the PDP-11/03 microcomputer to supply as a real-time processor synchronized 
to the MLS time-frame. Also, all project personnel were needed on the 
simulation studies (angle-and DME-receivers) . The status of the 
experimental system development, as it stands, is summarized in 
Table 5. 

The immediate objective with respect to an experimental system 
is the specification of the functional characteristics of a mini-or 
micro-computer suitable to the computational load of the receiver 
algorithm as it becomes firm. The work done to date is nearly machine 
independent, and the experience obtained will facilitate the drawing 
of a computer specification when appropriate. 
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PART TWO 


DME-RECEIVER STUDY 
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PRECEDING PAGE BLANK NOT FILMED 
SECTION VII 
DME SIGNAL MODEL 

The DME signal model is initiated at the RF stage o£ the receiver . 

At this point, the signal consists of a direct path signal, a 

reflected path signal, and added noise. It can be modeled as: 

X(t) = R(t) cos (tt :-„t) + pR(t-x) cos (m (t-x) + 3-,) + n(t) (7.1) • 

where R(t) = direct path envelope 

p = amplitude of reflected signal relative to the direct signal 

t = time delay of reflected signal 

3-^ = phase shift of carrier wave upon reflection 

n(t) = receiver noise assumed to be covariant stationary, Gaussian, 

' and bandpass with spectrum centered at the RF carrier, 

X(t) = R(t) cos (<u t) + p R(t-x) cos (co t + B 9 ) + n(-t) (7.2) 

c c 2 - 

- R(t) cos (m -t) + p R(t-x) rcos(g_) cos (a> t) - sin (3 0 ) sin 
c 2 c 2 

(w t)] + n(t) 
c 

where 3 0 = 3-, - w x 
2 1c 

Because of the assumed properties of the receiver noise process, 
it can be expanded into quadrature components with respect to the 
carrier frequency, oi^. 

n(t) =■ n c (t) cos(a> c t) - n_ (t) sin (to t) -(7.3) 

consequently X(t) can be written 

X(t) = [R(t) + p R(t-x) cos ( 3.0 + n (t)] cos (m t) - 

2 • -c c 

[p R(t-x) sin ($„) + n,-(t)] sin(w t) (7.4) 

2 s c- 

X(t) = X (t) cos(w t) + X (t) sin(m t) (7.5) 

c c s c 


53 



After heterodyning, the IF signal is 

Y^(t) = ^(t) cos (co Q t) + X g (t) sin(u o t) 

+ higher harmonics (7.6) 

With good tracking of the IF frequency, the output of the IF 
filter can be approximated as 

Y 2 (t) = Y c (t) cos(w o t) + Y g (t) sin(w o t) (7.7) 

where Y c (t) = X c (t) * h(t) 

Y g (t) = X g (t) * h(t) 

' h(t) is the impulse response of the lowpass analog of the IF filter 
*• denotes convolution 
It follows that the IF envelope is 

V (t) = CY c 2 (t) +Y s 2 (t)] 1/2 (7 . 8) 
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SECTION VIII 


DME SIMULATION 

The DME signal is simulated by generating the functions R(t), 

n (t) , and n(t) and combining them as in (7.4) to produce the 
c s 

quadrature components, X (t) and X (t). Each component is filtered 

c. s 

by separate but identical lowpass filters to produce the IF components, 
Y c (t) and Y g (t). The IF envelope is obtained by (7.8) and 
examined to determine when a threshold crossing has occurred. A large 
number of simulation runs are made (about 250) for - each set of 
multipath and noise conditions so that an approximate statistical 
average and distribution, function for the error will result. The 
error itself is then passed through a ten radian per second bandwidth 
filter to reduce random pulse— to— pulse errors. Another average 
and distribution function are obtained for the filtered error. 

8.1 IF Filters 

The quadrature components of the DME signal are filtered through 
a five-pole Butterworth lowpass filter. The effect of this filter 
is equivalent to that of the Butterworth bandpass filter used in the 
Hazel tine DME system. It is implemented in the simulation as a 
digital lowpass filter with a sampling frequency of 100 MHZ. Although 
a sampling frequency of this magnitude might be impractical in a real 
time situation, its use on the computer is justified since the 
response of the digital filter should be as close to that of the 
analog filter as possible. The bandwidth of the filter is 1.75 MHZ 
so that the simulation results will be directly comparable to the 
Hazeltine study [6]. 
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8.2 False Alarms 


In any type of receiver, there is a danger of mistaking an early 
threshold crossing due to noise as the' final crossing due to the 
DME pulse. .The rate of these "false alarms" varies with the’ 
proximity of the threshold level to the noise level and has the 
potential for causing a high negative bias in the error. Therefore, 
there is a ’need for some logic in the .receivers to recognize and 
eliminate some of these false alarms. 

.The primary factor that distinguishes a false alarm from the 
actual threshold crossing is the amount of time the signal stays 
above threshold. Samples of the- signal which are separated by a 
time constant (reciprocal of the bandwidth of the bandpass IF filter) 
are nearly uncorrelated so that the probability of the signal 
remaining above the threshold due to noise alone for a period of 
time greater than a time constant is very small. The half -amplitude 
pulse width is over ten times the length of a time constant so it 
should be possible to distinguish between false alarms and the 
actual threshold crossing due to the DME pulse. It is on this basis 
that false alarms are reduced in the simulation. 

8.3 Error 

Errors due to multipath and noise are combined here instead of 
being treated separately as in the Hazeltine and M.I.T. studies, [6], 
[7] . These two sources of error are not independent of each other 
since the signal to noise ratio (S/N) is affected whenever .multipath 
is present. Of course, if the threshold is placed at a very high point 
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on the pulse, multipath will be the primary source of error with noise 
errors being negligible.' The opposite is true if the threshold is 
placed at a very low point on- the signal, but since one tries to place 
the threshold at a point where there is an effective trade-off v 
between the two it is not realistic to analyze each separately. >,<: 
Previous studies [6] , [8] have assumed that the noise effects 
on the IF signal envelope, V(t) are Gaussian. The density function 
under these conditions is 

P = 1 exp [V (t) - V (t)] 2 / 2a 2 ] (8.1) 

V(t) (2Ha 2 )l/2 


where V Q (t) is the IF signal envelope uncorrupted by noise. 

This kind of assumption ignores the non-linear effects of the 
envelope detection process and as Rice [9] has shown,, the density function 
for the envelope is actually. 

M r ^tt2 


V(t) 


— ex P r(V 2 (t) + V o 2 (t))/2a 2 ] I o [V(t) V o (t)/a 2 ] ( 8 - 2 ) 


where I [V(t)V (t)/a 2 ] is the Bessel function, 
o o 

Rearranging terms, this density function is 


‘V(t) 


. V(t) 


exp j [(V(t)/cr) 2 + (V (t )/a) 2 ] I [(V(t)/a) (V (t)/a)] 


(8.3) 


If a(t) = V(t)/a and y(t) = V (t )/cr, then a new density function 
for the random process, a(t), will result. - 

Pa(t) = a(t) exp j [a 2 (t) + y 2 (t)] I q ta(t) y(t)} (8.4) 


This is a normalized density function for the envelope with the 
parameter y(t) representing the signal to noise ratio. 
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When (8*4) is evaluated at the threshold crossing time, a density \ < 
function for the amplitude variations in the I.F. envelope at this point 
in time results. Figure 14 shows this Rician density function for 
various values of y. The distribution is Rayleigh when y is zero 
and approaches a Gaussian distribution as y is increased, y in this 
case is the threshold to noise ratio (T/N) . 

A variation in the amplitude of the envelope can either shorten 
or lengthen the time it takes to reach threshold. If the error in 
the measurement of the arrival time of the pulse is considered to 
be the shift in the threshold crossing time from the ideal (the 
crossing time on an incoming pulse uncorrupted by interference from 
any form of multipath or noise) , then this density function also 
applies to the error. Positive and negative variations in the envelope 
amplitude cause early and late threshold crossings, respectively, 
so one would expect from the graph in Figure 14 that there would be 
a negative bias in the error at low T/N. Therefore, the assumption of 
a Gaussian distribution of error is only valid for relatively high 
T/N. 

8 . 5 Power Budgets 

The power budgets for the ground-based and airborne receivers 
are based on the landing pattern shown in Figure 15. The flight path 
is at a three degree angle with respect to a 14000 foot runway with 
the groundbased receiver at one end. Signal, attenuation is incurred 
along the flight path due to the antenna pattern [6] . This effect 
must be compensated for if accuracy is to be maintained as the plane 
approaches the runway. 
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For a cos-cos 2 pulse shape with a 3.5 p sec half -amplitude 
pulse width, the ratio of pulse peak to RMS amplitude is 3.7 decibels 
(DB) . The average power allowable for this pulse shape under the 
ICAO Annex 10 specification on effective radiated power in adjacent 
channels is 55 DB. This is obtained by integrating the pulse spectrum 
over all frequencies and comparing it with the power allowed in a 
0.5 MHZ band, 2 MHZ from the carrier. The peak is then 58.7 DBM. 

The assumed power budgets follow. The signal to noise ratios in 
the assumed budgets are defined as the ratio of peak signal level to 
RMS noise level and are defined at the input of the receiver before 
any power loss occurs (antenna losses are neglected) . Under these 
circumstances the noise power used is 4 kTB + NF instead 
of kTB + NF which is used to describe the available noise power in 
many cases. 


AIR TO GROUND POWER 


ERP (peak) . 58.7 

Transponder Antenna Gain 8.0 

Path Loss (18000 Ft) 107.0 

Peak Signal (18000 Ft) -43.3 

Noise (4kTB + NF) -87.5 

NF == 15 DB 

S/N (18000 Ft) 47.2 

S/N (16000 Ft) 44.2 

S/N (15000 Ft) 38.2 


BUDGET [6] 

DBM 

DB 

DB 

DBM 

DBM 

DB 

DB - 
DB 
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GROUND. TO AIR POWER BUDGET [6] 


ERP (peak) 

61.7 

DBM 

Path Loss 

107.0 

DB 

Peak Signal 

-45.3 

DBM 

Noise (4 ktB + NF) 
NF = 15 DB 

-87.5 

DBM 

S/N (18000 Ft) 

42.2 

DBM 

S/N (16000 Ft) 

39.2 

DBM 

S/N (15000 Ft-)" 

33.2 

DBM 
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SECTIONS IX - XII 


The following sections give simulation results for both the 
threshold and delay and compare receivers in the ground-hased transponder 
and the airborne interrogator. Each section is broken up into sub- 
sections %?hich describe how the receiver perforins as each parameter 
is varied in turn with the others .remaining constant. .This allows 
one to extrapolate as to the performance of the receiver under a 
wide range of conditions. 

An error summary is provided at the end of each section to 
provide a more detailed accounting of the performance of each 
receiver. This includes the performance under varying multipath 
conditions at the signal to noise ratio outlined in the power 
budget and also with a 6 DB drop in the signal to noise ratio 
as might occur with a specular reflection. The threshold level 
for this summary is chosen to be within a range of values where 
the receiver provides the best performance with respect to noise 
and multipath errors. The error in the summary has been filtered 
through a 10 rad/sec lowpass filter. This process has very little 
effect on the mean but reduces a by a factor of approximately 2.2. 
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SECTION IX 


SIMULATION RESULTS - TRANSPONDER 
FIXED THRESHOLD RECEIVER 

The threshold level is assumed to be set at a constant voltage 
level in the fixed threshold receiver. This means that the position 
of the threshold with respect to the pulse shape changes with the 
strength of the signal. A subtractive multipath with a relatively 
small delay would lower the signal to noise ratio (S/N) of the pulse. 
The resulting shift of the threshold to a higher point on the 
signal would cause a late threshold crossing and bias the error in 
the positive direction. An additive multipath signal would cause the 
opposite effect and negatively bias the error. 

This type of receiver is assumed to be limited to the ground- . 
based transponder so that the air to ground power budget applies here. 

9 .1 Threshold Levels 

The performance of the fixed threshold receiver is greatly 
dependent on the threshold setting. For a threshold to noise 
ratio (T/N) of 6 DB, the mean error is negatively biased for both 
additive and subtractive multipath signals (Figure 16). Early 
threshold crossings are due to the proximity of the threshold to 
the noise level and multipath has little effect under these conditions. 
A subtractive multipath signal, however, can provide a slight 
improvement since it causes a positive shift in the mean error 
bringing it closer to zero . 

The noise errors are diminished and a positive shift in the 
error results when the threshold is raised to higher levels. For 
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subtractive multipath, it increases almost linearly and 
becomes positive as the threshold level is increased, For additive 
multipath, the error reaches a minimum and shifts toward the negative 
again. The graph in Figure 16 can be separated into regions where 
noise is the primary cause of error (up to T/N = 9DB) and where multipath 
is the primary cause of error (T/N = 12DB and above) . The region between 
these two is where an effective tradeoff between the too sources of 
error occurs. 

9.2 Error Distribution 

The density function for the error has been described in section 
8-4. Figure 17 shows the error distribution for a threshold to 
noise ratio of 18 DB and a subtractive half-amplitude multipath. 

This is a nearly Gaussian distribution with a narrow spread about 
the mean. A very different distribution of error results when the 
threshold level is lowered to 6 DB (Figure 18). The mean is shifted 
into the negative region and the spread of error is no longer nearly 
symmetrical about this point. The error has now approached a 
Rayleigh form with the error spread over a large range below the 
mean and concentrated in small range above the mean. 

When the error is filtered through a 10 rad/sec lowpass 
filter, the spread is reduced by a factor of two to three but the 
original shape remains. This makes it difficult to obtain one 
expression which accurately expresses the spread of error about the 
mean. Using the standard deviation, cr, as a measure disquises the 
fact that the error is not symmetrical about the mean in all cases , 
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Error (Ft.) 

Figure 17 Error Distribution, Fixed Threshold Receiver, T/N - 18 DB 




but should provide a good approximation for thresholds in the range 

/ 

of interest. 

9 . 3 Changes in the Signal to Noise Ratio 

A change in S/N shifts the position of the threshold on’ the 
pulse and therefore changes the threshold crossing time. S/N is 
affected by changes in the signal strength due to specular ground 
reflections and receiver to receiver gain variations due tq temperature, 
aging and other factors. 

Ground reflections have very short differential path delays 
with respect to the direct path signal and this makes them recognizable 
only as changes in S/N at the receiver. Lateral reflections (multipath) 
have longer delays and are treated as a separate problem. 

Figure 19 shows the changes in the error bias for ± 6 DB 
changes in S/N xrfiile under the influence of a subtractive half- 
amplitude multipath signal. The 6 DB change can cause" a shift in the 
error bias of up to 33 feet. This shift generally increases as 
the threshold is raised. 

9.4 Multipath Effects 

Multipath in the fixed threshold receiver causes a shift in S/N 
which in turn causes error as explained above. The magnitude of the 
error is dependent upon the multipath parameters with error increasing 
as amplitude increases or as the differential path delay decreases. 

The error is most severe at short delay times as shown in Figure 20. 
It levels out to a relatively small error at about 150 nanoseconds (ns) . 
In the noise-free case, one would expect the error to level out at zero. 
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Figure 19 Mean Error Vs. S/N, Fixed Threshold Receiver 





When noise is added in, there is a negative bias in the error, the 
magnitude of which depends upon T/N. At T/N = 6DB, the noise bias is 
large so that short multipath delay times can actually cause an 
improvement in the error. When T/N is raised to 12 DB, the noise ' 
bias decreases and the error is near zero for multipath delays above 
150 ns. 
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0° 
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Multipath 
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a (Ft) 

E (Ft) 

a (Ft) 

E (Ft) 

a (Ft) 

E (Ft) 

a (Ft) 

Amplitude 

0.2 

(DB) 

46 

-4.4 

6.2 

-3.5 

6.4 

-2,6 

6.5 

-3.5 

6.4 

0.2 

40 

20.7 

8.1 

22.3 

8.4 

24.5 

8.8 

22.6 

8.5 

0.5 

46 

-5.5 

6.0 

-3.7 

6.3 

-1.0 

6.8 

-3.5 

6.4' 

0.5 

40 

18.3 

7.7 

22.1 

8.3 

28.1 

9.5 

22.9 

8.5 

0.8 

46' 

-6.8 

6.2 

-3.6 

6.3 

0.5 

7.6 

-3.4 

6.4 

0.8 

40 

16.3 

7.4 

21.4 

7.9 

32.8 

10.5 

23.2 

8.9 


Table 6 Summary of Errors - Transponder, Fixed Threshold Receiver 



SECTION X 


SIMULATION RESULTS - TRANSPONDER 
DELAY AND COMPARE RECEIVER 

A delay and compare receiver compares the IF envelope with a 
delayed an slightly amplified version o£ itself to determine the 
arrival time of the DME pulse, . A diagram of this type of receiver -is 

shown in Figure 21 . The threshold crossing time is given by the 
negative going zero crossing of the difference signal, d(t), and the 

threshold level is set by the delay and gain parameters, t and k. 

A form of automatic gain control is inherent in this type of receiver 
since the input signal is being compared with itself. 

Before the pulse arrives, the receiver will essentially be 
comparing noise signals which are highly correlated with each other 
due to the short time delay. The delayed signal will also be 
amplified so that there is a high probability that the difference 
signal will be below zero during this period.- Therefore, the. false 
alarm rate for this receiver is more significant than that of the 
fixed threshold receiver under similar conditions and must be 
reduced as outlined in section 8 . 2 . 

10.1 Threshold Crossing Point 

Figure 22 shows the mean error as a function of the threshold 
crossing point. The negative bias in the error due to noise occurs 
at higher levels on the pulse than it does in the fixed threshold 
receiver under the same circumstances . Assuming that the receiver 
delay, t, is fixed, it is necessary to increase the gain, k, to 
lower the threshold. This results in an increase in the noise level 


75 




Figure 21 Delay and. Compare Receiver Signals 






Figure 22 Mean Error Vs. 

Compare Receii 


old Crossing Point, Delay and 
= 46 DB 


in the delayed signal which in turn causes an increase in early 
threshold crossings and contributes to this effect. 

It is possible to separate the graph in Figure 22 into a region 
where noise is the primary source of error (up to T/N = 24 DB) and 
a region where multipath is the primary source of error (T/N = 29 DB . 
and above) . The region between these two is where the minimum 
over-all bias in the error due to multipath and noise occurs. 

Subtractive multipaths signals cause the most significant errors 
and also reach an overall minimum in this region, so these are studied 
in more detail . 

10.2 Error Distribution 

The density function for the error developed in section 8.4 
does not apply directly to the error in the delay and compare 
receiver. The error density function applies only to the-direct 
IF envelope signal, V(t), and not to the difference signal, d(t). 

However, the overall effect is much the same. Figure 23 shows 
the error distribution for a subtractive half -amplitude multipath 
signal at a T/N of 25 DB. This distribution is in the Rician 
form and is similar to that found at a 6 DB threshold level in the 
fixed threshold receiver. A threshold level of 29 DB gives an 
error distribution which is more nearly Gaussian (Figure 24) . 

The most, significant difference between the error distributions 
in the two receivers is the difference in the threshold levels where 
the near Gaussian distribution is acheived. This point is of interest 
because it indicates the signal level at which the receiver becomes relatively 
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Figure 24 Error Distribution, Delay and Compare Receiver, T/N = 29 DR 



insensitive to noise effects. It takes an increase in the 
threshold level of approximately -12 DB for the' delay and compare 
receiver to approach the same, level of insensitivity as the fixed 
threshold receiver. This is significant because the higher threshold 
levels mean more susceptibility to multipath errors. 

10 . 3 Changes in the Signal to Noise Ratio 

Positive or negative shifts in the signal to noise ratio may 
be caused by specular ground reflections or receiver to receiver gain 
variations. The effect of these shifts on the error bias is most 
pronounced in the case of a drop in gain (Figure 25) . The' inherent 
automatic gain control of the delay and compare receiver does not 
allow the threshold level to shift with respect to the pulse shape 
as S/N changes. The result is that the greatest penalty is incurred 
when S/N drops since this increases noise errors. A 6 DB S/N 
increase causes the error bias to improve slighty. 

The variations in the error bias decrease as the threshold 
level is raised due to a lessening of noise errors. 

10.4 Multipath Effects 

The M. I.T. study [7] has shown that the error caused by 
multipath alone in this receiver is small when the multipath delay 
is small, increases to a peak, and then diminishes steadily as the 
delay time is increased with all other factors remaining constant. 
This is also true when multipath and noise errors are considered 
together (Figure .-26). However, in this case when subtractive 
multipath is involved, the point where the peak error bias occurs 
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Multipath Delay (ns) 


Figure 26 Mean Error Vs. Multipath Delay, Delay and Compare Receiver 
S/N = 46 DB 




is dependent on the threshold level. The error bias at the 25 DB 
threshold level peaks at a multipath delay of approximately 50 ns. As 
the threshold is increased to 28 DB the peak error point shifts 
forward to about 100 ns. 

The key to understanding this 'is .the fact that subtractive 
multipath causes a change in the signal to noise ratio and at the 
point where the threshold is crossed there is a relatively larger 
change when the multipath delay is small. Low signal to noise 
ratios cause more noise errors (Figure ' 25) as do low threshold 
levels (Figure 26) so that the combined effect of both of these 
causes a shift in- She peak error. 
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5.1 
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46 
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Table 7 Summary of Errors - Transponder, Delay and Compare Receiver 
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SECTION XI 

SIMULATION RESULTS - INTERROGATOR 
ADAPTIVE THRESHOLD RECEIVER 

The adaptive threshold receiver takes the DME pulse at the IF 
stage and uses automatic gain control to normalize it. The threshold 
is set at a constant voltage level below the pulse peak and so 
does not shift its position on the pulse due to varying signal 
strengths as in the fixed threshold receiver. 

The effect of the AGC is such that the noise level is increased 
whenever a loss of signal strength occurs as in subtractive multipath 
conditions . This is in contrast to the delay and compare receiver 
which is self-AGC’d and thus does not change the noise levels. In 
cases of severe signal loss, the noise level could potentially 
be multiplied to the point that it approaches the threshold level. 

There should be a limit on the range of the AGC to prevent this . 

11.1 Threshold Crossings 

Under the influence of additive multipath, the error bias is 
positive and relatively constant through a large range of threshold 
levels (Figure 27). Under these conditions, noise errors are 
insignificant due to the decrease in the noise level caused by the 
AGC and as a result there is no negative shift in the bias at the 
lower threshold levels. 

The noise levels are increased by subtractive multipath 
causing a subsequent negative shift in the error bias at low 
threshold crossing points (up to T/N = 15.5 DB) . The error reaches a 
minimum point in the 15.5-20 DB region and again shifts toward the negative 
at higher threshold points due to the larger error penalties caused 
by the multipath. 
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Multipath and noise are not necessarily the primary causes 
of error for this receiver as they were for the delay and compare 
and fixed threshold receivers. Any loss of signal strengh caii 
potentially cause large errors due to the multiplication of the 
noise level by the AGC. This dictates the investigation of a 
somewhat larger range of threshold levels than in previous receivers . 

11.2 Changes in the Signal to Noise Ratio 

The shift in the error bias due to a 6 DB drop in the signal 
to noise ratio is about 24 feet at a threshold level of 15.5 DB 
(Figure ' 28). This a larger shift than any encountered in the 
two transponder receivers. As the threshold level is raised the 
shift decreases and finally reaches' a value of about 2 feet at the 
23 DB level. This suggests that the increased penalties in error 
bias incurred at the higher threshold levels may be offset by less 
sensitivity to signal degradation. 

11.3 Multipath Effects 

When the multipath delay is increased with all other factors 
remaining constant, the error increases to a peak at about 300 ns 
(Figure 29).- The adaptive threshold receiver is therefore sensitive 
to a larger range of multipath delays than either the fixed threshold 
or the delay and compare receiver. 

For multipath delays between 0 and 300 ns, the peak value of 
the pulse is approximately constant and as a result the AGC compensation 
is the same throughout this range. The effect on the leading edge 
of the pulse, however, is greatly dependent on the multipath delay 
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Figure 28 Mean Error Vs, S/N, Adaptive Threshold Receiver 
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and the AGC compensation. • The first portion of the leading edge 
is unaffected by the multipath so that when the AGC normalizes 
the pulse, it either sharpens or flattens this part of it. This 
causes the threshold crossing time to be pushed either backward 
or forward, respectively, from the ideal crossing time (Figure 30). 
This effect increases as multipath delay increases and the direction 
of the error shift is dependent only on the phase of the multipath. 
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tj - Ideal Threshold Crossing Time 

t 5Q - Crossing Time for 50 ns Multipath and AGC 

fc 100 - Crossing Time for 100 ns Multipath and AGC 


A - Direct Signal with No AGC 
B - Signal with 50 ns Multipath and AGC 
C - Signal with 100 ns Multipath and AGC 
D - Direct Signal with AGC 


Figure 30 Shift in Threshold Crossing Time with Multipath Delay, 
Adaptive Threshold Receiver 
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Table 8 Summary of Errors - Interrogator, Adaptive Threshold Receiver 



SECTION XII 


SIMULATION RESULTS - INTERROGATOR 
DELAY AND COMPARE RECEIVER 

I 

The delay and compare receiver used in the interrogator is 
identical to that used in the transponder. The conditions under which 
it operates are different, however, and are outlined in the ground-to- 
air power budget. 

12.1 Threshold Crossings 

The increased noise level in the interrogator results in a 
larger region where noise is the dominant source of error (Figure 31) . 
The minimum error bias with respect to both noise and multipath 
occurs i_n the 0.13 to 0.17 level range on a normalized pulse 
which translates into a range of threshold to noise ratios between 
22 and 25 DB. 

12.2 Changes in the Signal to Noise Ratio 

A 6 DB drop in the signal to noise ratio causes an increase 
in the error bias of 9 to 12 ft. depending on the threshold level 
(Figure 32). This degradation in performance is caused by an 
increase in noise errors. A 6 DB gain in the signal to noise ratio 
results in a slight improvement in each case. 

12.3 Multipath Effects 

The error caused by multipath in the delay and compare receivers 
under these conditions is larger than it is in the transponder. 

The higher threshold levels that are necessary here are the reason 
for this, making multipath and noise errors the major constituent in 
the combined error. 
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Figure 31 Mean Error Vs. Threshold Crossing Point, Delay and Compare 
Receiver, S/N = 40DB 
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Multipath errors peak at delay times of approximately 100 ns 
for threshold levels between 23 and 25 DB (Figure 33) . One would 
expect a shift of the peak point to shorter delay times and an 
increase in error at lower thresholds as in the delay and compare 
transponder. 
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SECTION XIII 


OVERALL STUDY CONCLUSIONS 

Angle-Receiver Study 

The integrated LOE/Kalman filter receiver algorithm tested in 
simulation as generally superior to the threshold receiver. Specifically, 
in the crossing-multipath scenario, primarily a test of tracking 
performance, improvement ratios (in peak absolute error) ranged to 
20:1 and better in certain situations involving high multipath 
interference. In the in-beam multipath and representative landing 
scenarios the optimal receiver superiority was confirmed, though less 
dramatically, partly because of the element of multi-path 
acquisition present in these runs. 

A distinct disadvantage of the optimal receiver is its complexity. 

The non-adaptive receiver (of the same structure), evaluated as a 
subop timal alternative, retained some of the superiority of the optimal 
receiver in multipath environments at a fraction of the computational 
demand. This suggests a carefully drawn compromise of performance 
and complexity might result in a computationally more efficient 
algorithm offering most of the principal benefits of the optimal 
receiver demonstrated. This problem area along with multipath 
acquisition . (identification) have been included in our plans for 
next year's effort. 

DME Study 

Under the assumed operating conditions of the transponder, the 
fixed threshold receiver seems to provide marginally better performance 
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than the delay and compare receiver. The fixed threshold receiver 
can have considerable immunity to both multipath and noise effects 
if the noise level is relatively low and the threshold is set at the 
proper point above this level. A disadvantage of this receiver is 
its sensitivity to changes in S/N which may be caused by specular 
reflections, receiver to receiver gain variations, and other causes. 

The delay and compare receiver has an inherent automatic gain control 
and is insensitive to these effects. 

The adaptive threshold receiver used in the interrogator performs 
poorly under any condition which reduces the input signal amplitude. 

The AGC under these conditions multiplies .the noise level and increases 
noise errors. The AGC also causes this receiver to be susceptable 
to multipaths with a large range of differential path delays. The 
delay and compare receiver with its inherent AGC can provide performance 
superior to the adaptive threshold receiver under all of these conditions. 
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ANGLE-RECEIVER INNOVATIONS STUDY 


In the scan data processor, new observations data are entered 
via a random process represented by the J-vector w(u|q) with representative 
element w_. , defined in (3.25), and repeated here with index j suppressed: 


A 

w = 


u 




I + qh 


where 


u = q + 2n c s/ q + n c 2 + n, 2 


q is a real number > 0 


(A. 1) 


(A. 2) 
(A. 3) 


n c ,n are independent Gaussian random variables with 

mean 0, variance 0.5 . (A. 4) 

The results of a simulation study of the first and second-order 
statistics of w (A.l) are given in Table A.l. The sample size was 
1000 points; the quantity RI in the table is on independent variable 
equivalent to twice the q parameter in (A.l) above. The autocorrelations 
shown are really values of the sample correlation coefficient, having 
been normalized to the appropriate^ sample mean square value. 

Conclusions dram are as follows: 

1. The sample mean (MEAN) is much less than the sample rms yalue 
•(WKMS) for all 5 RI values used and also it seems, as a 
random variable, to be well dispersed about zero; hence, it 
seemed plausible that 

<w|q )> = 0, independent of q (A. 5) 

and this, conclusion ‘was ''dr awn. 

* , " . „ < ‘ 3 
* \ »** * ‘ 

2. The sample correlation coefficients for non-zero shifts are 

• much less than unity for all 5 RI values used, suggesting that 
a sequence of w-values with q fixed is a white process; the 
whiteness property was assumed to extend to the more general 
non-fixed q case. 
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The further observations, concerning the sample mean square value 
(WMS) , that 

1. WSS Z 1, for RI = 0 

2. WMS 1 , for RI large 

RI 

suggested a tenative approximation formula, in terms of q, as follows: 

<^w 2 |q^ = h(q) : ~ 2 q (A. 6) 

The results of a more extensive simulation, involving 10,000 samples 

RHO 

and values of q ~~ 2 ~) f rom 0 to 50 (corresponding to S/N = 34db), 
are shown in Table A. 2 and in Figures A.l and A. 2, comparing plots 
of the sample -mean square value and the approximation (A. 6). The 
error in the approximation peaks at about 20% for q = 2 (RHO — 4) 
and seems in an average sense to be asymptotic to zero for smaller and 
larger values of q. The approximation (A. 6) was employed in the 
scan data processor with good results. 
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ABSTRACT 


This paper addresses an estimation problem in which a landing air- 
craft uses ground-transmitted microwave information to determine its azi- 
muth angular position 0 ( t ) relative, to a fixed reference. State estima- 
tion is used to lower the mean-square error in estimates of 0(t) produced 
by an envelope processor in the airborne receiver', 0(t) is modeled as 
part of the state of a linear dynamic system driven by white Gaussian 
noise of unknown covariance. The envelope processor estimates become 
linear observations of the state corrupted by additive Gaussian noise of 
known covariance. Adaptive. Kalman filtering is examined as a means of 
computing estimates of 0(t) having minimum mean-square error. Adaptive 
filtering methods are found in the technical literature which work for 

N 

systems where che noise is stationary. They are then modified for use in 
the aircraft position estimation problem, where the noise statistics are 
time varying. The adapti-ve filters are tested in a digital computer sim- 
ulation, where e{t) is updated according to aircraft motion along an 
unknown flightpath. Several of the adaptive filters work very well, 
though not significantly better than suboptimal estimators of less 
com pi ex’ 1 ’ ty * 
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CHAPTER I 


INTRODUCTION 

This paper describes the application of state estimation theory to 
an aircraft landing problem where the system model is incompletely 
defined. In general the problem requires estimation of the state of a 
linear dynamic system driven by white Gaussian noise with unknown covari- 
ance. The state is observed by a linear function of the state corrupted 
by additive white Gaussian noise. When all model parameters are known, 
the optimal minimum variance estimator becomes the Kalman filter [1, pp. 
228-229], [2, pp. 195-201]. However, when the model noise covariances 
are unknown, the optimal estimator cannot be achieved, and some subopti- 
mal approach must be employed. Several adaptive Kalman filtering methods 
from the literature are examined in this paper as possible solutions to 
the aircraft landing problem. 

Before giving a formal description of the state estimation problem, 
let us first provide a background description of the aircraft landing 
problem. A more rigorous problem definition can then be presented, along 
with a proposed course of solution. 

Background on Microwave Landing System 

The problem examined in this paper is part of an airborne receiver 
study for the Microwave Landing System (MLS). The MLS, developed by the 
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Federal Aviation Administration (FAA) and the National Aeronautical and 
* ' / 

Space Administration (NASA), provides electronic guidance in an air ter- 
minal area for landing aircraft [3], [4]. The system enables an 
approaching aircraft to compute its position in space relative' to a fixed 
ground reference. The required coordinate information is derived by the 
aircraft's receiver from ground-transmitted microwave signals. 

Let us establish a cartesian coordinate system, with its origin at 
the stop end of the runway. Referring to Fig. 1-7. A, the runway center- 
line forms the X axis, while the Z axis is normal to the ground plane. 

We also establish a spherical reference system centered at the same 
origin.. At time t, the aircraft's position shall be defined by the fol- 
lowing spherical coordinates: 

r(t) = direct path distance from the origin to the 
aircraft. 

s(t) = azimuth angle from the X axis to the projection 
onto the ground plane of a ray from the origin 
to the aircraft. 

<j>(c) = elevation angle from the ground plane up to this 
ray. 

The MLS enables the aircraft to compute these three coordinates. We 
restrict ourselves in this paper to considering only the azimuth angle 
0(t). First we present a brief description of the azimuth channel in 
the MLS. 

An antenna, located at the coordinate origin, electronically • scans 
a ±60° azimuth coverage sector with a narrow fan-shaped microwave beam. 
The beam is narrow only in the azimuth sense (1° between -3 db points). 




B. Scanning Beam Boresight vs. Time 



o ^ 0®) 

C. A Noiseless Signal Envelope for 0 (t ) — 20 


Figure I - 1 
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while being wide in elevation coverage. The beam is scanned in a 
"TO-FRO" fashion, as shown in Fig. I-l.B: the beam boresight e^(t) 

starts at +60° azimuth, moves at constant rate to -60°, holds there for a 
brief "dead time," and then moves back. This scanning procedure takes 
12.2 milliseconds and repeats at a 13 1/3 Hz update rate [3, pp. 1-10, 

11, 27]. 

We define t^ as the time at the start of the kth scanning proce- 
dure. We also assume that e(t) is constant at 0 (t^) during the 12.2 mil- 
lisecond duration of the scan. During this time the scanning beam 
signal, viewed at the input to a receiver at the aircraft, is amplitude- 
modulated, having large amplitude when the boresight 0^(t) is near e(t^). 
The envelope-detected signal, shown in Fig. I-l.C, has two -pulses: one 

which peaks wrren s^(t} = 0(t^) during the "TO" scan; and the second which 
peaks when & A (t) = 0(tj,} during the “FRO" scan. As seen in Figs. I-l.B, 
C, the time differential between the centroids of the t wo pu lses is 
directly related to the value of e(t k ). The aircraft can therefore 
determine its azimuth angle by receiving and envelope detecting the 
ground-transmitted signal and measuring this time differential. 

This scheme for computing ©(t^) runs into difficulty when we real- 
ize that the received signal is corrupted by front-end noise in the air- 
borne receiver. This front-end noise produces random distortions in the 
envelope so that any attempt to estimate © ( ) from envelope information 
will have random errors as well. 

An optimal envelope processor has been designed which, given the 
noise-corrupted signal envelope for the kth scanning interval, computes 
an estimate of s(t k ) which minimizes mean square error. The 
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envelope-detected IF signal is sampled in the vicinity of the two large 
pulses during the 12.2 millisecond scanning process. Then, during the 
"down time," before the next scanning process begins, the envelope sam- 
ples are sent to a minicomputer. Here a "locally optimum estimation" 

t 

algorithm computes an estimate of [5, pp. 8-29], [6,pp. 4-28„ 

52-61]. Using a stochastic model for the signal envelope, this algorithm 
provides an estimate of 0(’t^), given the envelope samples for the kth 
scanning process, which is optimal in terms of minimum mean square error. 

The optimal envelope processor estimates e( . based only upon the 
envelope samples taken on the kth scanning interval. At a 13 1/3 Hz azi- 
muth update rate, we would expect the effects of thermal noise upon the 
signal envelope to be independent between consecutive scans. The error 
in consecutive estimates should therefore be independent as well. On the 
other hand, the true angle 0(t, ) cannot change appreciably between scans 
for a large aircraft. In plotting a time sequence of estimates we 
therefore expect to see random fluctuations- about a slowly changing 
mean. 

Since the estimates change much more rapidly than the true azimuth 
angle, it seems reasonable that the estimate of 0^) could be improved 
by averaging it with past estimate values. This would produce a new 
estimate based on all past envelope information and not just that 
obtained on scan k. This is the objective of the work presented in this 
paper. Adaptive Kalman filtering is examined as a means of producing an 
estimate of ©(t^) having a smaller mean square error than that of the 
envelope processor estimate. 
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Thesis Overview 

The above-stated problem is presented mathematically in Chapter II, 
where a stochastic system model is derived. The Kalman filter requires a 
state variable model, where an nth-order linear dynamic system is driven 
by white noise of known covariance.' We therefore model ©(t^) as part 
of the state of such a system. , The discrete Kalman filter is presented 
in Chapter III as the optimal estimator of ©(t^), given the state model of 
Chapter II. It is also shown that the Kalman filter requires knowledge 
of the plant noise covariance, which is unknown in our problem model. 
Adaptive Kalman filtering is therefore studied as a suboptimal estima- 
tion approach in which the unknown Kalman gain is estimated from measure- 
ment information.. 

Several candidate adaptive-filtering methods from the literature 
are presented in Chapter IV. Each filtering scheme is developed under 
the assumption of stationary noise. In Chapter V we modify each of the 
candidate filters to work for our specific problem, where the unknown 
noise covariances are time varying. 

The adaptive Kalman filters are tested in digital computer simula- 
tion in Chapter VI. This testing proceeds in two stages. First the 

r 

assumed stochastic model of Chapter II is simulated .with additively cor- 
rupted measurements of the state sent to a candidate adaptive filter. 

The error in the filter’s estimate e(t^) is plotted as a function of 
time. If a candidate filter performs well here, it is then tested in a 
second simulation phase where the stochastic model assumption is 
removed. ©(t^) 1S now updated deterministically as the aircraft moves 
along a prescribed flight path. Additively corrupted measurements of 
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0(t k ) are again sent to the candidate filter which then computes the 

■•N 

estimate 0(t^,). That adaptive filter is sought which minimizes the mean 
square error in ©(t^,). A conclusion is given in Chapter VII. 

A’; 



CHAPTER II 


PROBLEM DEFINITION AND MODEL DEVELOPMENT 

In this chapter we offer a more rigorous problem description and 
then develop a stochastic model describing the evolution of the azimuth 
angle sCt^). This model is then used in subsequent chapters to develop 
an estimate of e{tjJ, based on all past envelope information. 

Problem Definition 

Before formally describing the estimation problem, we place some 
mild restrictions on the aircraft's azimuth coordinate and its estimate 
produced by the envelope processor. We first change notation, using O(k) 
instead of to represent the azimuth coordinate at the start of the 

kth scannina interval. 

w / 
Let us assume that the aircraft is making a landing approach along 

some prescribed flight path unknown to us. As the aircraft moves along 

this path, let its azimuth angle be given by 

e(k) = f(k.) •• (Il-l) 

While we do not know the relation f(*) s we shall assume that it is a 
member of a known "class" of functions representing evolutions in e(k) 
for typical landing approaches. For example, we can limit the aircraft's 
maximum air speed or minimum radius of turn. More is said about this in 
Chapter VI. ' 


8 
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We define y(k) as the estimate of ©(k) produced by the envelope 
processor using the locally' optimum estimation algorithm. This estimate 
is unbiased and can therefore be given by 

y(k) = 0(k) + v(k) C I I- 2 ) 

- 

where v{k) is a zero-mean additive error term with covariance R(k). R(k) 
is computed by the locally optimum estimation algorithm so that we know 
its value [6, p. 5]. The probability distribution of v(k) is unknown. 
Here., we assume that it is Gaussian. This does not appear to be an 
unreasonable assumption insofar as we would intuitively expect the error 
density to be symmetric about a single.mode at zero. ‘Also, the Gaussian 
assumption makes the state estimation problem to follow mathematically 
tractable. We therefore write the probability density of v(k) as 

p[v(k)] = N[0, R{k)j ( I I — 3 ) 

We hereafter denote the first-order density of an n-dimensional Gaussian 
process x(k) with mean m(k) and covariance P(k) as 

p[x(k)] = iN[m(k), P(k)] (II-4) 

where 

_n_ 

N[m(k},P(k)3 ? ( 2ir ) 2 |P(k)| 2 exp{-l[x(k)-m(k)] T P' 1 (k)[x(k)-m(k)]} (II-5) 

If x(k) is white or uncorrelated in time, we write: 


p[x(.k)] = WN[m(k) , P{k)] 


(II-6) 
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The envelope processor error v(k) is produced by the effects of front-end 
noise on the IF signal. These effects are independent from one scan 
interval to the next, so that v(k) is uncorrelated. We therefore write 

p[v(k)] = WN[0, R(k}] ( I I - 7 ) 

As mentioned in Chapter I, the envelope processor estimate y(k) 
uses only the envelope information from the kth scan interval. Our • 
objective is to develop an estimate of e(k) based on Y^, the set of all 
past values of y(k) : 

Y k 4 {y(D, y(2), y(k)} (n-8) 

For a slowly changing azimuth angle there is a high correlation between 
©(k) and efk +• 1), while the estimates y(k) have uncorrelated errors from 
scan to scan. As mentioned in Chapter I, we intuitively expect to 
improve the estimate y(k) by averaging it in some way with past values. 
This could be viewed. as a low-pass filtering approach. 

Let us consider a stochastic state model, driven by noise, as a 
representation for the evolution of 0(k). Given a valid state model, we 
could then, by treating the estimates y(k) as observations of the state, 
produce a new state estimate which minimizes error in some mean-square 
sense. This is the approach taken here. 

We now offer a formal problem description. Given in the problem is 
an unknown azimuth coordinate ©(k) described by (II-l), where f(») is a 
member of a known class of functions. Also given is the set Y^ of past 
envelope processor estimates. The estimate y(k) has zero-mean, white 
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Gaussian error with known covariance, described by (I 1-2) and ( I 1-7 ) - 
The objective can be stated as follows: using an assumed stochastic 

model for the evolution of 0(k), develop an adaptive Kalman filter which 
estimates e(k) so as to minimize mean-square error in the estimate. Sev- 
eral adaptive Kalman filtering methods are obtained from the literature 
and modified for use in this problem. Each candidate filter is tested in 
computer simulation, with error in the estimate of 0 ( k) being the quant- 
ity of interest. 

Stochastic Model 

The Kalman filter requires a state variable model where an nth- 
order 1 i'near system is driven by white Gaussian noise. Vie therefore use 
such a. stochastic model in representing the evolution of the azimuth 
angle e(k). In order to keep the resulting Kalman filter computationally 
feasible we elect to use a two-dimensional model where the angular accel- 
eration is set equal to white .Gaussian noise: 



p[u(t)] = WN[0, S(t)] (11-10) 

Our decision to model acceleration as white noise provides us with the 
lowest order stochastic model for which both e(t) and ©(t) can be esti- 
mated. As shown in Chapter III, e(t) is used to linearly extrapolate the 
estimate of ©(t) between measurement times. 

In using the noise process u(t) to model ©(t), we must relate the 
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noise covariance S(t) to the actual value which 0(t) takes on. Since 
S(t) is the expected value of the square of u(t), we set it equal to the 
square of the actual acceleration: 

S(t) = e 2 (t) (II-T1) 

Of course, in estimating G(t) we do not know o(t), since the only availa- 
ble information is the set of envelope processor estimates. S(t) is 
therefore unknown in our model, at least from the aircraft's point of 
view. 

As stated in the problem definition, we are interested in estimat- 
ing e(k), the value of e(t)at the start of the kth scanning interval. 

We therefore obtain a discrete-time representation of the state model in 
{II-9}. Ler us first replace ( I I -9 } with a more general state equation: 

x(t) = Ax(t) + Gu(t) (11-12) 

where x(t) is a general state vector driven by a vector noise process 
u(t). (11-10) can still be used to describe u ( t) . A general discrete 

state model is given by: 

x(k) = $(At)x(k - 1) + Fw(k - 1) (11-13) 

p[u)(k - 1)] » WN[0, Q(k - 1)] (11-14) 

where At - t^ - t^ _ ^ . (11-14) becomes the discrete equivalent of 

( 1 1 -9 ) with x(k) representing x(t)L ^ f when we use the following 

Z \ 

transformations [7, pp. 60-61 , 72-75] : 
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$(At) = exp{AAt} (11-15. A) 

t 

rQ(k - l)r T = / k 4(t k - t)es(t)G T j T (t k - t)'dt (11-15. B) 

t k - V 

Substituting the A matrix of ( I 1—9) into (11-15, A), we obtain for our 
model : 

#(At) = 

We assume At to be constant, so that $(At) also becomes a constant and is 
written hereafter as $. For evaluation of (II-15.B) we assume that S(t) 
is constant at S(t^) over the limits of integration. This seems reasona- 
ble, as O(t) cannot change appreciably during one scan period at a 
13 1/3 Hz update rate. Moving the scalar S(t^) outside the integral and 
changing variables we obtain: 

T At T T 

rQ(k - l)r - S(t. ) / <Kt)GGV ( T )dT (11-17) 

K 0 ■ 

We then substitute for G and $(t) from (II -9) and (11-16) and evaluate 
the integral: 

(11-18) 


rQ(k - 1 )r T = s(t k ) 


1 

¥ t 


lat 2 

At 


1 At 
0 1 


' (11-16) 


The linear at term in the matrix dominates for small at, as is the case 



14 

for our problem, where At = .075 (the period of a 13 1/3 Hz update rate). 
We approximate the other terms as zero: 

rQ(k - l)r T = 


0 0 

o Atsct^) 


(11-19) 


With this approximation Q(k - 1) becomes a scalar, so that the state 
equation becomes 


efk) 


1 At 


0(k - 

1) 


0 


= 







0(k) 


0 1 


0(k - 

1) 


1 


p[a(k - 1)] = WN[0, Q(k - 1)] 


(11-20) 
( 11 - 21 ) 


Q(k - 1) = AtS(t k ) (11-22) 

Since S(t) is unknown in the continuous-time model, Q(k - 1) is unknown 
as well. 

We can see from (11-20) that the effect of our assumption in 
(11-19) is to add noise only to the velocity e(k), so that e(k) becomes 
piecewise linear between measurement updates. If we were to assume a 
constant acceleration between times t^ _ -j and t^, e(k) would be 
described by : . 

e(k) = o(k - 1) + Ate(k - 1) + ^At 2 Q(k - 1) (11-23) 

1 2 " 

In (11-20) we have discarded the nonlinear term ^At e(k - 1), assuming it 
to be negligible. For reasonable landing approaches we do not expect 0 
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to exceed 0.1°/sec^ [6, p. 40]. For a 13 1/3 Hz update rate, At is 75 
milliseconds, and the error in neglecting the nonlinear term is always 
less than ^(.075 sec)^(.l°/sec^) , or 2.8 x 10"^ degrees. Since 0.01° is 
given as a desired r.m.s. error specification, (11-20) is a valid 
approximation. 

We now have a discrete time state model describing the evolution of 
e(k). We recall from ( 1 1 -2 ) and ( 1 1 -7 ) that the envelope processor esti- 
mate y(k) equals Q(k) plus a Gaussian error term. We can therefore view 
y(k) as a linear observation, or "measurement" of the state corrupted'by 
additive noise : 

+ v(k) > (11-24) 

Finding the value of 9(k) now becomes a state estimation problem. We 
must estimate- the state of a linear dynamic system excited by white 
Gaussian noise of unknown covariance, given linear measurements of the 
state corrupted by additive white Gaussian noise of known covariance. 


y00 = D 0] 


©00 

©00 



CHAPTER III 


THE DISCRETE KALMAN FILTER 


In this chapter we examine the discrete Kalman filter, which is the 
optimal estimator for the assumed state and measurement models in our 
problem. As previously mentioned, the optimal estimator cannot be used 
here, as the state noise covariance is unknown. The optimal estimator is 
of use, however, in obtaining the suboptimal solutions to follow, and 
provides a lower bound on error performance. 

Let us first provide the state and measurement equations in concise 
form. From (12-7), (11-20), (11-21), and (11-24) we have: 


e(k) 


1 At 


e(k 

- 1) 


— 

0 


= 





+ 


e{k) 


0 1 


e(k 

- 1) 


1 


1) 


(HI-1) 


y(k) - [ i 


0 ] 


o(k) 

0(k) 


+ v(k) 


(III-2) 


p[<o(k - 1)] = WN[0, Q(k - 1)] ( 1 1 1 — 3 ) 

p[v(k)] = WN[0, R(k)] . ( 1 1 1 -4)* 

where Q(k - 1) is unknown. Equations ( 1 1 1 -1 ) and ( 1 1 1 -2 ) describe a spe- 
cific member within a general class of linear systems given by: 


16 
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x(k) = *x(k - 1) + rw(k - 1) ( 1 1 1 - 5 ) 

y(k) = Hx(k) + v ( k ) (III-6) 

where x ( k ) is an D-dimensional state vector, y(k) is an m-dimensional 
measurement of x(k), and where to(k - 1) and v(k) are Gaussian noise vec- 
tors of dimension r and m. The noise sequences for the general case are 
still described by (III —3) and (III-4), where Q(k - 1) and R(k) are now 
symmetric, non-negative definite matrices of respective dimensions r x r 
and m x m. 

Let us consider the general system of (111-5) and (III-6). We 
assume in this chapter that the state noise covariance Q(k - 1) is known. 
Our objective is to estimate the state x(k) given the set of all past 
measurements: 

Y k A {yO). y(2), • • • y(k)} (m-7) 

Let us define x(k[j) as an estimate of x(k), given Y.. We are concerned 

J 

with finding the optimum state estimate x(kjk). 

In order to have a quantitative measure for optimality we define a 
performance index, or loss function J(kjk): 

J(k|k) = E{ [x ( k) - x(kjk)] T W[x(k) - x(kjk)]} ( I I 1-8) 

where W is an n x n symmetric, non-negative definite matrix. When W is 
diagonal J(kjk) becomes a weighted sum of the mean-square errors in the 

A 

elements of x(k|k). We define the optimal state estimate of x(k) as that 
estimate x(kjk) which minimizes J(kjk). It can be shown that J ( k | k ) is a 
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member of a class of loss functions which are minimized by that estimate 
given by the conditional mean of the state given all past measurements; 

, >*s 

i.e., the optimal estimate x(k|k) becomes 

x(k|k) = E{x(k)|Y k } C I I I— 9 ) 

[1 , pp. 227, 231], [2, pp. 147-148], This is true for any non-negative 
definite W. 

We again reference the general linear system and 

still assume that the noise covariances are known. It is well known that 
the optimal estimate of x(k) for this system, the conditional mean which 
minimizes J(kjk), is given by the discrete Kalman filter [1, pp, 228- 
229], [2, pp. 195-201], The Kalman filter 'is described by the following 


equations [2, p. 201]: 

x(k jk - 1) = $x(k - 1 [k - 1) ' (III-TO) 

P(k{ k - 1) = $P(k - 1 jk - T)* T + rQ(k - l)r T (III-ll) 

K(k) = P(kjk - 1 )H T [HP(k|k - 1)H T + R(k)]" 1 (III-12) 

x(k|k) = x(k|k - 1) + K(k)[y(k) - Hx(kjk -1)] (III-13) 

P(kjk) = P(k|k - 1) - K(k)HP(k|k - 1) (III-14) 


where x(k[k - 1) denotes the optimal predicted, or extrapolated estimate 
of x(k) given Y^ while x(kjk) is the optimal updated estimate, using 
all measurements up to the present time. The term P(kjj) is the error 
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covariance in the estimate x(kfj): 

P(kjj) = E { [x ( k ) - x ( k | j ) ] [x ( k ) - x(k|j)] T } (III-15) 

K(k) is the Kalman gain, which determines the weighting given the present 
measurement y(k) in computing x(kjk). Note that K(k) is not a function 
of measurement values, so that x(k)k) is a linear estimate of x(k)'. 

Returning to our original problem of ( 1 1 1 - 1 ) to (III-4), we define 
the state estimate x(k|k) by 

x(k|k) = [S(kjk), e(k|k)] T (III-16) 

We seek. to minimize mean-square error in 0(k]k), so that our performance 
index becomes 

Pl(kjk) = E { [ 0 ( k ) - e(k[k)] 2 } (III-17) 

From (III-S) we see that PI ( k [ k ) is a special case of the general loss 
function J(k|k) where W is the diagonal matrix diag{l, 0}. Therefore the 
Kalman filter produces the optimal estimate for our problem when the 
noise covariances are known. We now give the Kalman filter equations for 
our specific model: 

©( k i k - 1) = 0 ( k - 1 jk - 1) + Ate(k - Ijk - 1) (HI-18) 

P ]1 (kjk-1 ) = P-j-j (k-1 ]k-l ) + 2AtP 12 (k-l [ k-1 ) + At 2 P 22 (k-l { k-1 ) (III-19) 

P ]2 (kjk - 1) = P 12 (k - Ijk - 1) + AtP 22 (k - 1 Jk - 1) (III-20) 
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P 22 (k|k - 1) * P 22 ( k - 1 Ik - 1) + Q(k - 1) 

K-j ( k ) = P 11 (klk - l)/CP 11 (k!k - D + R(k)3 
K 2 (k) - P 12 (kjk - 1 )/[P-j -j ( k 1 k - 1) + R(k)l 

e( k { k) » e(k|k - 1} + K 1 (k)[y{k) - e(k|k - 1)] 

£(k|k) = e(k - l|k - 1) + K 2 (k)[y(k) - e(k|k - 1)1 
' P 11 (klk) = P-,-, (k|k - l)[l - K-j (k)3 
P 12 (k[k) = P 12 (k{ k - 1)D - K-j(k)] 

‘ P„( ktk) = P 22 (k(k - 1) - K 2 (k)P 12 (kjk - 1) 


(III-21 ) 
( 111 - 22 ) 
(111-23) 

(II 1—24) 

$111-25) 
(III-26) 
(III— 27) 
(111-28) 


• Le t us consider the Kalman filter for the general case of an n- 
dimensional state vector with scalar noise' and measurements. From (III- 
12) we can see that the optimal gain matrix K(k) becomes an n-dimensional 
vector. Some of the properties of this gain vector are useful for the 
adaptive Kalman filter development of Chapter IV. 

We first note that for stationary noise where Q(k - D and R(k) are 
scalar constants, the gain K(k) as well as the error covariances 
P(k|k - 1) and P( k j k ) reach constant steady-state values. If Q(k - D 
and R(k) are constant or slowly varying, K(k) is a function of the ratio 
Q(k - 1 )/R(k) . Furthermore, K(k) can be specified if only one of its 
elements are known; all remaining entries are deterministic functions of 
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the one known element. The dependence of the gain on Q and R as well as 
the functional relationships between gain elements can be observed numer- 
ically, but we cannot usually find closed-form expressions for such pro- 
perties, especially when n is large [8, p. 274]. Another useful property 
of the gain for the scalar measurement case is given by Alspach [8, p. 
272]: 


0 < HK(k) < 1 


(III-29) 


We offer a proof for (III-29), but first we rewrite the general 
Kalman gain equation (III-12): 

W(k) = HP ( k [ k - 1 )hT + R ( k ) (II 1—30 ) 

K(k) = P(k[k - 1 )H T W _1 ( k) (III-31) 

P ( k. | k - 1) and R(k) are non-negative definite matrices, so that W(k) is 
non-negatiye definite as well, having the same form as a covariance 
matrix. More is said about this property of W(k) in Chapter IV. Alspach 
notes that W(k) can be written: 

W(k) = [I - HK(k)] _1 R(k) (III-32) 

(8, p. 270]. We can prove this by starting with the right-hand side and 
substituting: 


[I - HK(k)]' 1 R(k) 

= [I - HP ( k | k - l)H T ir 1 (kj] _1 R(k) 


{II 1—32 . A ) 


(III-32.B) 
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= W ( k ) [W ( k ) - HP( k [ k - l)H T ] -1 R(k) 

(II 1-32 - C ) 

- W(k)[R(k)] _1 R(k) 

(II 1—32 . D) 

= W(k) 

(II 1-32 . E) 


(III -32 . B ) is obtained by substituting for K(k} from (III -31 ) , while 
(II 1-32 - D) is obtained from (II 1-30) . For the scalar noise and measure- 
ment case R(k) and W{k) become positive scalars, and (III -32) can be 
rewritten: 

1 - HK(k) = R(k)/W(k) (III -33 ) 

Since R(k) and W(k) are positive, their ratio cannot be negative. From 
(II 1—30) we know that W(k) must be greater than or equal to R(k), so that 
the ratio in (III-33) cannot exceed unity. We' therefore have: 

0 < 1 - HK(k) < 1 (III-34) 

or 

0 < HK(k) < 1 { III-35) 

For the two-dimensional aircraft problem we can show that the gain 
of the filter given by (lII-18)-(III-28) is a function of Q ( k - 1)/R(k). 
We can also obtain a closed-form expression for the second gain element 
in terms of the first. We first need to express the error covariances 
P(k|k - 1) and P(kjk) in terms of K(k), From ( I 11-22) we can write: 

„ P-j-j ( k|k - 1) = k 7 (k)R(k)/[l - K-j (k)3 


(III-36) 
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Substituting this expression for P^(k[k - 1) into (II 1-26 ) we obtain: 

P n (k[k) * K 1 (k)R(k) (II 1-37 ) 

Substituting (III -36 ) for P^-j(kjk - 1) into (III-23), we have: 

Sr 

„ 

P-j 2 (k [ k - 1} = K 2 (k)R(k)/[l - K-j (k)] (HI-38) 

And substituting this expression for P-j 2 ( k | k - 1) into ( 1 1 1 -27) we have: 

P-, 2 ( k I k ) = K 2 (k)R(k) (III— 39) 

We now make the assumption that the noise covariances are slowly 
varying in time. This seems reasonable, as Q(k - 1) is determined by the 
acceleration eit^} and R(k) depends on the signal -to-noise ratio. Nei- 
ther of these quantities can change appreciably between scan intervals at 
a 13 1/3 Hz' update rate. We therefore approximate the error covariance 
P ( k J k - 1) as having the same value for two consecutive time periods: 

P ( k + 1 }k) : P(k]k - 1) (III -40 ) 

Using this approximation; we substitute (II 1-38) and (II 1-39) into (HI- 
20) and obtain: 

P 22 (kjk) = K 1 (k)K 2 (k)R(k)/[At(l - ICj (k))] (II 1-41 ) 

We can now represent the second entry of the optimal gain in terms of the 
first. Using the approximation of (111-40)", we substitute (III -36 ) , 
(III-37), (III-39), and (III -41 ) for the needed covariances into (III-19) 
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and simplify, obtaining: 


K 2 (k) = K 2 (k)/[At(2 - K-j(k))] (III- 42 ) 


We observe that K^(k) is a monotone-increasing function of K-|(k). Using 
{111-35} and (II 1-42 ) and noting that here HK(k) equals K^(k), we estab- 
lish bounds on the gain: 


0 < K-j (k) <1 
0 < K 2 (k) < 1/At 


(III —43 ) 


Keeping the assumption of (III-40), we find from ( I I I -21 ) and (II 1-28 ) 
that we have two expressions for the difference P 22 (kjk - 1).- P^Ckjk). 
Equating these we have: 

Q(k - 1)'= K 2 (k)P 12 (kjk - 1) (II 1-44 ) 

Substituting (III-38) for P-j 2 ( k | k - 1), we have: 

Q(k - 1 )/R(k - 1) = K|(k)/Cl - K-j ( k ) ] (111-45} 

And by using (111-42} for ( k ) we finally have: 

Q(k - 1)/R(k) = K^(k}/[At(2 - K 1 (k)) 2 (l - K-j { k } ) ] ‘ (111-45} 

Q(k - 1)/R(k) is clearly a monotone-increasing function of K^(k), ranging 
from zero when K-j(k) is zero to infinity when K^(k} is unity. Since 
Q(k - 1 }/R ( k) and K^(k) are both positive, we can infer that K^(k) is a 
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monotone-increasing function of Q(k - 1)/R(k) as well. 

We can make some intuitive observations from this about the optimal 
filter of {III-18]-(III~28) . From (III-18) we see that ©(kjk - 1) is a 
linear extrapolation of e(k - l|k - 1), based upon the velocity estimate 
e(k - Ijk - 1). We' then accept a new measurement y(k) of 0(k)-, and use 
the weighted difference between y(k) and ©(kjk - 1) to update our esti- 
mate to ©(kjk) in {XI 1-24) . The weighting applied to this difference is 
K^{k), varying from 0 to 1 . We have just found K^(k),to be a monotone- 
increasing function of Q(k - 1)/R(k). We might view Q(k - 1)/R(k) as the 
ratio of uncertainty in our knowledge of the state x(k) to uncertainty in 
the measurement y(k). When this ratio is low, indicating high confidence 
in our estimate of the state, K(k) is small, so that y(k) has little 
effect on the new estimate x ( k j k ) . When this ratio is high, we have 
greater confidence in our new measurement y(k). K(k) becomes large, and 

A 

y(k) has more weighting in determining x{k|k). Of course when 
Q(k - 1}/R(k) approaches infinity, we have no prior knowledge of the 
state: x(k - Ijk - 1) gives no information about x(k). K-j(k) becomes 

unity, causing e(kjk) to become y(k). 

Before leaving the subject of optimal filters, let us study the 
effects of using a Kalman filter with incorrect or suboptimal gain. This 
would be the case if incorrect values were used for the noise covariances 
Q(k - 1) and R(k). Assume that the filter of (III-10) and ( II I -1 3) is 
implemented, with a suboptimal gain K(k). From (III -5 ) , (III-10), and 
(II 1-13) we can write the error in the state estimate: 
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x(k) - x( k | k) = $x(k - 1) + ro)(k - 1) - §x(k - l|k - 1) 
- K(k)[y(k) - H$x(k - 1 [k - 1}] 


(II 1-47) 


Substituting for y(k) from ( I I I -6) we have: 


x{k) - x( k J k) = $[x(k - 1) - x(k - Ijk - 1)] + r<u(k - 1) 

- K(k)[H$x(k - 1) + Hrw(k - 1) * v(k) - H$x{k -'Ijk - I}] 

(II 1-48) 


x ( k ) - x ( k | k ) = [I - K(k)HMx(k - 1) - x(k - Ijk - I)] 

(III-49) 

+ [I - K(k)H]r w (k - 1) - K(k)v (k) 


We recall that w (k - 1) and v ( k ) are samples of white sequences and are 
independent of each other. Since x(kjj) is a combination of measurements 
through y(j), we can make the following assertions: 

■ k < i 

[x(k) - x(k| j)] is independent of w(i): " (III-50) 

j < i 


all k 

[x(k) - x(k|j)j is independent of v(i): (I I 1-51 } 

j < i 


Therefore all three terms in (II 1-49) are independent, and we write the 
suboptimal error covariance: • . 


P(kjk) = E{[x(k) - x(k j k)][x(k) - x(k|k)] T } 


(III- 52 -A) 



= [I-K(k)H]$E{[x(k-l )-x(k-l | k-1 )][x(k-l )-x(k-l [k-1-)] V[I-K(k)H] T 


+ [I - K(k]H]rE{o)(k - l)J(k - l)}r T [I - K(k)Hj T + K(k)E{v(k)v T (k)}K T (k) 


(II 1-52 . B ) 


Simplifying this expression, we obtain: 


P { k } k ) = [I - K(-k)H][>P(k - Ijk - 1)* T + rQ(k - l)r T ][I - K(k)H] T 

(III-53) 


+ K(k)R(k}K T (k) 


Equation ( I 11—53) has been iterated until steady-state is reached 
for our problem of (III-l )- ( III— 4) with stationary noise. Figure III-l 
shows P-j-j(k|k), the steady-state mean-square error in e(kjk), as a func- 
tion of suboptimal gain. P-^(kJk) versus K-j(kjk) is plotted for three 
ratios of Q(k - 1)/R(k). Note that for each case P^(k|k) is minimum for 

a 

the optimal gain and then rises to R(k) as fCj(k) approaches unity. We 
can see that as long as the suboptimal gain is near or above the optimal 

A 

value, the estimate e(k|k) will have a lower mean-square error than y(k). 
When the suboptimal gain becomes less than the optimal value, however, 
P-|^(kjk) rises rapidly, approaching infinity as the gain goes to zero. 
Here the suboptimal filter diverges. The gain is so small that insuffi- 
cient weighting is given to the most recent measurement y(k) in updating 
the state estimate x(k|k). Too much emphasis is placed on old measure- 
ment information, so that the filter cannot follow the true dynamics of 
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CHAPTER IV 


ADAPTIVE KALMAN FILTERING: THE STATIONARY NOISE CASE • 

In Chapter III we introduced the discrete Kalman fil terras the 
optimal state estimator for the linear system described by: 


x(k) = <t>x(k - 1 ) + rw(k - 1 ) 

(iv-U 

y(k) = Hx(k) + v(k) 

(IV— 2) 

p[u(k - 1)] = WN[0, Q(k - 1)] 

(IV— 3) 

p[v(k)] = WN[Q, R(k)] 

' (IV-4) 


We made the assumption, however, that the noise covariances were known. 
Reviewing the Kalman filter equations (III-l 0)— (III-l 4) we note that both 
Q(k - 1) and R(k) are needed for computing the optimal gain K(k) and the 
error covariances P(k[k - 1) and P(kjk). If either Q{k - 1 ) or R(k) is 
unknown, as is the case for our problem, the gain K(k) cannot be found. 

We could implement the Kalman filter equations by substituting estimates 
for the unknown noise covariances, but the performance of the resulting 
estimator could be highly suboptimal if these estimates are poor. 

In this chapter we introduce adaptive Kalman filter as a suboptimal 
estimation scheme when the noise covariances are unknown. The adaptive 
Kalman filter takes the form of the optimal filter: 
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x(k[k - 1) = sx(k - 1 } k - 1 ) ( IV-5) 

x ( k ( k ) = x(k|k - 1) + K(k)[y(k) - Hx(kjk - 1)] (IV-6) 

where K{ k) is an estimate of the unknown optimal gain K(k). The subopti- 
mal K(k) is a function of the measurements Y^: we use the measurements 

to either estimate K(k) directly or to estimate the unknown noise covari- 
ances for use in the Kalman equations (III-lO)-(III-l 4) . The adaptive 
filter is therefore a nonlinear estimator, unlike the optimal Kalman 
filter, which is linear since the gain K{ k) is independent of the 
measurements. 

Here we present three adaptive Kalman filtering methods for the . 
stationary noise case from the literature, as well as a simpler intuitive 
scheme. Each method assumes the system model of (IV-1 ) — ( IV— 4 ) , with both 
Q ( k - 1) and R( k) unknown and constant. Some methods assume a model of 
general dimension, while others assume scalar noise and measurements . In 
this chapter we present the development of each method for the general 
stationary case. In Chapter V we modify the adaptive filters to work 
when the noise covariances are time-varying and apply them to our speci- 
fic problem of (III— 1 ) to ( I II— 4) . 

The Innovations Sequence 

Before presenting a development of the adaptive Kalman filtering 
methods, let us first define the innovations sequence. This concept is 
useful in the filtering developments to follow. 

We first recall the general Kalman filter equations (1 11-10)- (HI- 
14). Specifically, the updated estimate is given by (III-13) ; 
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x { k } k) = x{k[k - 1) + K(k)[y(k) - Hx(kjk - 1)] ■ (IV-7) 

We shall, define the innovations residual v(k) by: 

' 3 . 

v(k) = y(k) - Hx(k|k - 1) ' . (IV-8) 

The time sequence of these residuals is called the innovations sequence. 
We can show that the residual v(k) is actually the error in the optimal 
predicted estimate of the measurement y(k) given _ -j . From Chapter 
III we recall that the optimal, least-mean-square error estimate of y ( k) 
is the conditional mean of y(k) given the available measurements. From 
(II 1-9) we write the optimal predicted estimate as: 

y(k|k - 1} = E{y(k)jY k _ } } 

= E{Hx(k) + v{k)|Y k _ -j} 

■ HE{x(k)jY k ; ^ + E{v(k)} (IV-9.C) 

Equation (IV-9.B) results from substituting { I V-2 ) for y(k) } while the 
second term in (IV-9.C) results from noting that v{k) is from a white 
sequence and thus independent of Y k _ ^ . Recalling that the optimal pre- 
dicted state estimate is given by the Kalman filter, and that v(k) is 
zero-mean, we have: 


(IV-9.A) 

(IV-9.B) 


y(k|k - 1) = Hx(k|k - 1) 


(IV-10) 


where x(k|k - 1) is from the Kalman equation (III— 10). The innovations 
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residual v(k) is therefore the difference between the measurement y(.k) 
and its optimal predicted estimate. This error is then multiplied by the 

A 

Kalman gain K(k) and used to correct x ( k J k - 1) in ( IV— 7 ) to produce the 
optimal state estimate x(kjk). 

The residual v(k) can easily be shown to ’be zero-mean: 

E{v(k)} « E{y(k) - Hx(kjk - 1)} (IV-ll.A) 

= HE{x(k) - x(k| k -r 1)} + E{v(k)} = 0 (IV-ll.B) 

(IV-ll.B) is obtained by substituting (IV-2) for y(k)J It equals zero 
because the Kalman estimate x(k|k - 1) is by definition unbiased, while 
v ( k } is zero-mean. We can also find the innovations covariance W(k): 

W(k) = E{v(k)v T (k)} (IV-12.A) 

= E{[H(x(k) - x { k | k - 1)) + v(k)][H(x(k) - x(k|k - 1)} + v(k)] T } 

(IV-1 2.B) 

A 

From ( I I 1-51 ) , [x(k) - x ( k ( k - 1)] and v(k) are independent, so that (IV- 
12.B) becomes: 

W(k) = HE{[x(k) - x ( k f k-1 )][x(k) - x(kjk-l + E{v(k)v"**(k)} 

(IV-i 3 .A) 

= HP ( k f k - 1)H T + R(k) (IV-1 3} 

From (III-12) we see that the optimal filter computes W(k) in 
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finding the Kalman gain K(k). The gain was written as a function of W(k) 
in (II 1-30 ) , (III -31 ) , although no physical interpretation was given for 
W(k) at the time. We sometimes find it convenient to express the Kalman 
gain and updated state estimate of (III-12), (TII-13). in terms of the 


innovations sequence: 

v(k) = y(k) - Hx ( k i k - 1) (IV-14) 

W(k) = HP(kjk - 1)H T + R(k) (IV-15) 

K(k) = P(kjk - 1 }hW) (IV-16) 

x(k|k) = x ( k f k - 1) + K(k)v(k) (IV-17) 


The innovations residual v(k) can be shown to be Gaussian. Assum- 
ing the initial state x(0) to be Gaussian, we can see from (IV-1) that 
x(k) is Gaussian, since it is a linear combination of Gaussian random 
variables. Similarly from (IV-2) we see that y(k) is the linear combina- 
tion of Gaussian random variables and must be Gaussian as well. Finally 
we recall that for the optimal filter x( k | k - 1) is a linear combination 
of the measurement values in _ -j , and is therefore Gaussian. Since 
v(k) is a linear combination of y(k) and x(kjk - 1), it must also be 
Gaussian : 

p[v{k)] = N[0, W(k)J . CIV-18) 

We have already established in (IV-8)-(IV-10) that o(k) is the 
zero-mean difference between y(k) and its optimal estimate y(kjk - 1). 
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Since y(k|k - 1) is the conditional mean E{y( k) J ^}, we can represent 
the conditional density p[y(k)|Y^ _ as the density of v(k) with its 
mean displaced to y ( k | k - 1): 

p[y(k)jY k _ ^ = N[Hx(k[k - 1), W(k)] (IV-19) 

The innovations sequence becomes important when we realize that it 
is an obtainable measure of the estimator's performance. .From (IV-13) we 
know that the innovations covariance W(k) is directly related to the pre- 
dicted estimate error covariance P(k|k - 1). The derivation of (IV-13) 
makes no assumptions of filter optimality, so that this relation holds 
whether, the filter gain is' optimal or not. (Of course P(k|k - 1) and 
P(k[k) are computed for the optimal case, and W(k) is not needed). For 
our specific aircraft model of ( 1 1 1-1 ) — (III-4) , equation (IV-13) is given 
by (III- 22 ) : 

W(k) = P 11 (k[k - 1) + R(k) (IV-20) 

For a constant R(k), we see that the innovations covariance rises and 

a 

falls with the mean-square error in ©( k j k - 1). W(k) should therefore be 
minimum when the state estimator is optimal. 

We can solve for W(k) as a function of suboptimal gain K(k) for our 
specific system model of (III-l ) - ( 1 1 1 -4) when the noise is stationary. 

We have already obtained the error covariance P(kjk) as a function of 
suboptimal gain K and constant covariances Q and R in Chapter III by 
iterating (II 1-53) until steady-state is reached. From (IV-20) we have 
W(k) as a function of P ( k ( k - 1). We can therefore express W(k) in terms 
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of K by relating P{k|k - 1 ) to P(k|k) for the suboptimal filter. Using 
(IV-1) and (IV-5) we write the error in the predicted estimate: 

x(k) - x ( k f k - 1) - *[x(k - 1) - x(k - 1 Jk - 1)] + r<o(k - 1) (IV-21) 

Recalling from (III-50) that [x(k - 1) - x(k - Ijk - ))] and oj(k - 1) are 
independent, we have: 


P(k jk - 1) = E{[x(k) - x(kj k - l)][x(k) - x(kfk - 1)] T } .(IV-22.A) 

. = i>p(k - 1 [k - Us 1 + rq(k - l)r T _• (IV-22) 

This is the same relation as (III-10) for the optimal filter, which is 
given by (III-l 5)— ( II I-2T ) for our specific aircraft problem. To obtain 
W(k) we only need P,^(lc|k - 1), given by (III-19): 

P 11 fkjk-7) = P 11 (k-1 [k-T) + 2AtP 12 {k-l |k-1) + At 2 P 22 (k-l Jk-1 ) (IV-23) 

We cari therefore find W(k) as a function of K, Q, and R .by first obtain- 
ing P(k|k) from the steady-state solution of ( I 11-53) and then applying 
(IV-21) and (IV-20). 

Figure IV-1 shows plots of W(k) versus suboptimal gain K-j with con- 
stant Q(k -:1) and R(k) for aircraft system model (III-l ) - ( 1 1 1-4) . We 
note that W(k) is minimum when K-j equals the optimal gain K-j , This is 
expected, since W(k) is the sum of P,^(kjk - 1) plus R(k), and both 
P^(k|k - 1) and P^(kjk) are minimized when the filter gain is optimal. 
We also note that W(k) rises toward infinity as approaches zero; here 
P^(kjk - 1) and P^(kjk) are 'both approaching infinity as the suboptimal 




«> 


37 


filter diverges. We observe this same effect in Figure III-1. 

.The properties of the innovations sequence presented here are of 
immense value in the adaptive Kalman filtering developments which follow. 
It has been shown that the innovations sequence contains all new state 
information obtained by the measurements [9, p. 176],. In, addition to 
the properties stated above, Mehra shows that the innovations sequence is 
white for the optimal filter and correlated when the filter gain becomes 
suboptimal [9, p. 177], This property is not used by the adaptive fil- 
tering methods presented here. 

We are now ready to present methods of adaptive Kalman filtering 
for general stationary noise problems. 

The Method of Sage and Husa 

Let us assume the general system model of (IV-1 ) - (IV-4) , where the 
state vector x ( k ) has dimension n and the measurement y(k} has dimension 
m: 


x(k) = $x(k - 1 ) + rid(k - l ) 


(IV-24) 

y{k) = Hx(k) + v(k) 


(IV-25) 

P U(k - 1)] = WN[0, Q(k - 1)] 

- 

(IV-26) 

p[v(k)] = WN[0, R(k}] 


{ IV-27) 


The noise covariance matrices Q(k - 1) and R(k) are unknown constants, 
and shall be written hereafter as Q and R. We recall from (II-8) that 
is the set of all past measurements , and we define Xj, as the set of all 
past state values: 
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Y k = (y(l), y(2) , * • • y(k)} (iv-28) 

X k = (x(l), x(2), * • . x(k).} ( I V— 29) 


Sage and 'Husa develop maximum a posteriori (MAP) estimates x(k[k), 
'Q(kJk), and R(kjk) which maximize the conditional probability density of 
the unknowns given the measurements . They actually address the more gen- 
eral problem where the noise terms of ( IV— 26} and (IV-27) have unknown 
means to be estimated as well. This more general procedure i.s not appli- 
cable for our problem, however, and is not covered here. The reader is 
referred to the works of Sage and Husa for a description of their general 
method [iO], £11 J. 

Let us form the a posteriori density of the unknown states and 
noise covariances given the measurements, i.e.., the conditional probabil- 
ity density of X k , Q,. and R given Y k : 


r v 

pla 


k 3 


Q, RJYJ = 


P[Y k |X k , 


Q, R]p[X, 
P[YJ 


Q, R] 


(IV-30) 


The right-hand side of (IV-30) is obtained from Bayes Law, where 
p[X k , Q, Rj is the a priori density of the unknowns given no measurement 
information. For maximum a posteriori estimation we need to find those 
values of x(k), Q, and R which maximize (IV-30), Noting that the denomi- 
nator p[Y k ] is unaffected by the choice of these values, we seek to maxi- 
mize the function: 


J(k) = p[Y k |X k , Q, R]p[X k , Q, R] 


(IV-31 ) 
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Sage and Husa assume that the unknowns are independent in the absence of 
measurement information, yielding: 

J(k) = p[Y k |X k , Q, RMX k ]p[Q]p[R] ' (IV— 32) 


They next assume that the a priori densities of Q and R-are *uni form 
between some known limits, 'For example: 


P[Qij]~i- 0 . _ n > 

I n jMAX Hi jMIN 

L o 


°ijMIN - Q ij- Q ijMAX 
otherwise 


( I V— 33) 


All we know about theijth element of Q is that it lies somewhere between 
Qi JMIN an< ^ ^ijMAX" ^ values between these limits are equally likely. 

Of course, if we have no information on how large Q and R become, we can 
allow the limits to approach infinity (Q-j must positive for diago- 
nal elements, as Q is non-negative definite). As long as we remain 
within the allowable limits on Q..and R. . the densities p[Q] and p[R] are 

ij ij 

constant and do not affect the maximization of J(k). We therefore write: 


J(k) = Cp[Y k |X k , Q, R]p[X k ] (IV-34) 


where C is a constant. 

We now solve for J(k) in terms of its component densities. p[X k ] 
can be expressed: 

p[X k ]'= p[x(k), X k . ,] (IV-35.A) 


* P[x(k)jX k _ 1 ]p[X k _ -,] 


(IV-35) 
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where the last step results from the definition of conditional probabil- 
ity. We then substitute (IV-24) to write: 

P[ x (k) | _ -j] = p[$x(k - 1) + Tco(k - l)|X k _ ^ (IV-36) 

Given X k _ ^ , x(k - 1) is known. No new information about rto(k - 1} is 
obtained, since from (III-50) to(k - 1) and x(j) are independent for 
j < k -.1. p[x(k)jX k _ -j] therefore assumes the density of rto ( k - 1) shifted 
in mean by §x(k - 1 ) : 

p[x(k)jX k _ -j] = N[$x(k - 1), rQr 1 "] (IV-37) 

Reapplying (IV-35) and the substituting (IV-37) we have: 

■ k 

PCX. ] = p[x{0)] n p [x ( j ) j X . ,] • (IV-38.A) 

K : j=l J • 

k T 

= N[x(0), P(O)] ir N{>x(j - 1), rQr 1 ] ' (iv-38) 

j=l 

We now reapply the definition of conditional probability to obtain: 

p[Y k |X k , Q> R] = P[y(k) s Y k . _ 1 |X k> Q, R] (IV-39.A) ' 

- p[y(k) 1 Y k _ r X k , Q, R]p[Y k _ 1 ]X k , Q, R] . (IV-39) 

Using (IV-25) for y(k) we write; 

p[y{k)jY k _ ]> V Q» R] “ P[Hx(k) + v(k}|Y k _ r X k , Q, R] (IV-40) 
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Given we know the value of x(k). Knowledge of the conditioning vari- 
ables gives no new information about v(k), which is independent of X k and 
y(j) for j < k. (IV-40) therefore becomes the density of v(k) with mean 
displaced by Hx(k) : 


p[y(k)[Y k _ v X k , Q, R] = N[Hx(k) , R] - ''' - (IV-41) 


In the same manner used to obtain (IV-38) we reapply (IV-39) and then 
substitute (IV-41 } : • . 


p[Y k ] x k , Q, 


k 

R] = n 
. 0=1 


p[y(J) 


y j 


- i 


V Q, 


R] 


(IV-42) 


k 

= n N[Hx(j) , R] (IV-43) 

0=1 


We can now solve for J(k) by substituting (IV-38) and (IV-43) into 
(IV-34): 


k k 

J(k) = CN[x(o), p(o)] n N[#x(j - i), r Qr T ] r N[Hx(j), r] 

0=1 0=1 

(IV-44.A) 


• = CN[x(0),P(0)] 

.k JL . . -1 ' ‘ • -.1 .. 1 

x -n ( 2 tt ) 2 1 rQr T J 2 exp{-kx(jHx(j-l )] T (rQr T ) [x(j)-*x(j-l)]} 

j=l ■ c 


; ^ . .lii I 1 . • • - ’ 

x . n (2n) 2 1 R { 2 exp{-i[y(j)-Hx(j)] T R -1 [y(j)-Hx(j)]> 
0=1 


(IV-44) 
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Incorporating into C all components of J ( k ) which are unaffected by 
Q, and R and therefore constant for the maximization procedure, we have: 


d(k) = CjrQr T | 



■— k . ' 

2 exp t [x{j)-$x(j-l)] T (rQr T ) 


-1 ■ . . 

[x(j)-«x(j-l)] 

• ( IV-45) 


- 1 s Cy(a )-H x(j)I] t r"' 1 Cy( j }-Hx( j)] T > 
^j=i 


We now have J(k) as a function of X^, Q, and R. Needed at time k 

A A A 

are the values x(kjk), Q(k|k), and R(kjk) which maximize J(k).. Sage and 
Husa solve this problem by using a "discrete maximum principle" [10, p. 
770], Here we offer an alternate approach' yiel ding the same results. 

We first maximize J(k) with respect to Q and R. 

To obtain the- MAP estimate Q(k[k) we rewrite- (IV-45): 


k _i - - 

J n (k) = c|rqr T | 2 ex P (4 z [x(j)-$x(j-l)] T (rQr T ) [x(j)-*x(j-i]} • (IV-46) 


2 j=l 


where C contains all factors of J(k) which are not functions of Q. We 
choose to redefine J in terms of ln(d) and maximize this function 
instead. This is allowable, since 1 n ( - ) is a monotone increasing 
function: 


k -1 

Jq(k) = -kl n J rQr T ( - s [x( j ) - $x(j-l )] T (rQr T ) [x(d) - $x(j-i)] 


(IV-47) 


Note that Jq(k) is a function of rQr T , so that we cannot estimate Q 

T 

directly. Substituting P for rQr and z(j) for [x(j) - $x(j - 1)], we 
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have: 


J(P) = -klnjP| - z z T (j)P _1 z(j) (IV-48) 

3=1 

We seek that value Pq of P which maximizes (IV-48]. 

Let us define e as a scalar arbitrarily close to zero, such that 
Pq + eA represents a small deviation in P from Pq. Since P represents 

T 

fQr , both Pq and A are non-negative definite. We can write: 

. J(P 0 + eA) = J(P 0 ) - 6J(P 0 , e) (IV-49) 


Obvious.lv, 5J = 0 when e =' 0. Since J(P) is maximum at P Q , SJ cannot go 
negative and is minimum at P = Pq, or at s = 0: 


~{3J(P 0 , e) >L=0 = 0 


( I V— 50 ) 


We now must obtain a functional relation for 6 J(Pq, s) in order to 
solve for Pq. Using (IV-48), we write: 


J(P 


0 + eA) = -klnjpQ + eA | - S Z T (j)(P 0 + eA)" 1 z(j) 

j i 


(IV-51) 


where 


|P„ + eA| = |P„(I + eP n _, A) | - |P„||I + cP/’aI (IV-52) 


We use the approximation, valid for small e: 


1 + eB ; 1 + e trace(B) 


( I V- 53 ) 
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to write: 


| P 0 + eA l = \ P qH ] * e trace(P Q _1 A)] 
We next use the matrix identity [12, p. 79]: 

(B + cr 1 = B -1 - + C" 1 ) B" 1 

to write: 


{ IV— 54) 


( IV— 55 ) 


-1 

(P Q + eA)" 1 = Pq" 1 - P 0 ^[P 0 _1 + (eA)' 1 ] Pg" 1 (IV-56) 

For small e- (eA) - ^ » P^ , yielding the approximation: 

■ (P Q + eA)" 1 Z P 0 _1 - ^Pq'^Pq' 1 ( IV— 57 ) 

Now we substitute (IV-54) and (IV-57) into (IV-51): 

d(P a + sA) = -kln|P 0 | - kln[T + e trace(P 0 _1 A)] 

k k (IV “ 58) 

. - I z T (j)P n ' 1 z(j) + e Z z T (j}P n _1 AP '‘‘zU) 
j=l U j=l U 

We recognize the. first and third terms of (IV-58) as J C Pq ) from (IV-48). 

The second and fourth terms give us 5 J(Pq> e), so that we have: 

-1 k 

^-{6J(P 0 ,e)}=-k[l+e trace(P 0 “ 1 A)] [trace(P 0 -1A)]+_l z T (j)P Q " 1 AP 0 " 1 


(IV-59) 
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From (IV-50), we let this become zero for e = 0 , yielding: 


k trace(P 0 _ 1 A) = _s^z T (j)P ( 3 ~ 1 AP 0 ' 1 z( j) 


(IV-60) 


Using the matrix identity for symmetric, non-negative definite 8 [1, p. 
231]: 


x T Bx = trace(Bxx T ) 


(IV-61) 


we have: 


k tr 5 ce(P G “^A) = I trace[P 0 ” 1 AP 0 “ 1 z(j)z T (j)] 


( I V- 62 ) 


trace(kR 0 _ 1 A) = traceCP^ 1 APq -1 z( j )z T ( j )] 


( I V- 63 ) 


Equating trace arguments and simplifying: 


P n = l I z(j)z T (j) (IV-64) 

U k j=l 

Recalling from before (IV-48) that z(j) = [x(j) - 4 >x( j - 1)] and that P Q 

J 

is that value of rQr that maximizes Jq(k), we have: 

Ic 

rQ(k| k)r T -72 [x(j) - <j>x(j - l)][x(j) - <px(j - 1)] T (IV-65) 

\pl 

The state x(j) is unknown and will be replaced by the optimal smoothed 
estimate of x(j) given the measurements Y^: 
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rq( k J k)r T = TT z [x(jjk) - $x(j - l|k)][x(j|k) - $x(j - Ilk)] 1 ( IV-66) 
K j=l ' 

We obtain the MAP estimate R(k|k) in similar fashion, defining 
U R (k) those factors of J(k) in (IV-45) which are functions of R: 

JL k 

J R (k) = C | R j 2 exp{4 z [y(j) « HxUJjV 1 [y(j) - Hx(j)]} (IV-67) 

T 

J R (k) has the same form as Jg(k) in (IV-46), with R replacing rQr and 
[y ( J ) - Hx(j}j replacing [x(j) - ax(j - 1)]. Maximization .of J R (k) with 
respect to R should therefore yield an estimate R(kjk) of the same form 
as rQ(kjk)r' in (IV-65) : 

- R { k ; k ) = ~z [y( j } - Hx(j)][y(j) - Hx(j)] 1 (IV-68) 

We agai-rr replace x(j) with the optimal smoothed estimate x(j[k): 

~ 1 k ^ -r 

R ( k ]' k ) = -r- Z [y ( j ) - Hx( j |k)][y{j) - Hx(j|k)] (IV-69) 

K j=l 

We make the assumption that the optimal MAP estimates of (IV-66) 
and (IV-69) are very nearly equal to the true noise covariances: 

rQ( k | k)r T : rqr T , R(kjk) ; R (IV-70) 

A 

Under this assumption we can obtain the optimal MAP state estimate x(k|k) 

T 

from the linear Kalman filter of ( III-l 0)-(III-14) , where rQr and R are 
replaced by their MAP estimates. This assumption also allows us to 
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obtain the estimates x(jjk), needed by rQ(k | k)r T and R(k|k), from optimal 
1 inear smoothi ng. 

In using the Kalman equations (III-l 0)-(lII-14) to compute x(kjk), 
we must substitute R(k - l]k - 1) for R. x(kjk) requires R in the gain 

* 4 

equation (III-12), and since R(kjk) requires x{kjk) in (.IV-69), it does 
not yet exist. R(k - l|k - 1) is the best available estimate of R for 
computing x(k|k), and is therefore redefined: 


R(kjk-l) = R(k-1 jk-1 ) 


k-1 

j Z [y(j) - Hx(jjk)][y(j) - Hx( j|k)] T 


(I-V- 71) 


Sage, and Husa develop an estimation algorithm which uses (IV-56) 

- T 

and (IV-71) to compute FQ(kjk)r and R ( k j k - 1) for use in the Kalman 
filter. However the result quickly becomes complicated and impractical, 
because- of the 'need to process smoothed estimates x(j[k) [10, p. 762]. 
They then derive from (IY-66) and (IV-71) equations for computing subop- 

- T 

timal estimates rCj (kjk)r and R (kjk). These equations require only the 

A A 

estimates x(j}j - 1) and x(j|j) produced by the Kalman filter. We now 
present a development of their suboptimal mechod. 

In the suboptimal design to follow we first assume that the esti- 
mates for rQr T and R will be good enough that the Kalman filter using 
them will be nearly optimal. We therefore assume that the error covari- 
ances P(kjk - 1) and P(k|k) are computed correctly by the Kalman filter. 
We first replace R ( k ( k - 1) in (IV-71) with a suboptimal estimate 

A 

A ( k | k - 1) which does not require smoothing: 
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i k-i T 

. A(k|k - 1) = j- Up s [y(j ) - Hx( j ] j - 1 )][y(j) - Hx( j ] j - 1}] T 

(IV-72.A) 

i k-1 

= 2 v(j)v'(j) . ( IV— 72 ) 

3-1 

where v(j') is the innovations residua] given by (IV-8). To be unbiased, 

/V 

A ( k | k - 1} must have R as its expected value: - 

- 1 T 

E { A { k J k - 1)} * -j— -y Z E{v(j)v (j)} (IV-73) 

K " 3=1 

£'{v{ j)v T (j) } = W = HF( j j j - 1 )H T + R - , (IV-74) 

where W is- the steady-state innovations covariance, given by (IV-13). 
E{A(k|k - 1)} thus equals W, so that A(k|k - 1) is biase d. We note how- 

T 

ever that KP ( j j j - IjH 1 is computed in the Kalman filter's gain equation 
(III-12) , and can therefore be subtracted out of the summation term to 
produce an unbiased estimate: 

R (kjk - 1) = p-UpZ v(jK(j) - HP(j j j - 1)H* (IV-75) 

" 3=1 

We now replace rQ(k]k)r^ in (IV-66) with a new suboptimal estimate: 

~ i k „ A ~ „ j 

A(k|k) =~ Z [x(jjj) - $x(j)j - 1 ) ] [x ( j* 1 3 ) - $x{jjj - !)] (IV-76.A) 

J ^ 
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1 k T 

= 77 E [K(j)v(j)][K(j)v(j)] T (IV-76) 

where (IV-76) follows from noting that K(j)v(j) is the difference between 
estimates x(j|j) and x(jjj - 1) in (III-13). To find whether A ( k [ k) is 
unbiased we compute: 


E{A(k|k - 1)} = 1 Z K( j)E{v( j )v T ( j ) }K T ( j ) 

K j=l 

(IV-77) 

K(j)E{v(j)v T (j)}K T (j) = K(j)W(j)K T (j) 

(IV-78. A) 

- K(j)W(j)[P(jjj - l)H T W _1 (j)] T 

(IV-78. B) 

= K(j)HP(j|j - 1) 

(IV-78. C) 

= P(j|j - 1) - P(jjj) 

(IV-78. D) 

» ®P(j - 1 | j - 1)« T + rqr T - P(j|j) ' 

(IV-78) 


(IV-78. B) results from. substituting (III-12) for K(k), while (IV-78.C) 
results from noting that P(j|j - 1) arid W(j) are symmetric by definition. 
We then obtain (IV-78.D) from (111-14) and (IV-78) from (HI-11). A(k | k) 
is biased, but we note that $P(j - l|j - 1)® and P(jjj) are computed by 
the Kalman filter and. can be removed from the summation. We therefore 

T 

obtain the suboptimal estimate for rQr : 

rQ (k|k)r T * J- z K(jMj) J(j)K T (j) + P(j|j) - <&P(j - l|j - 1)* T 

S J=1 


( I V- 79) 
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We now write (I.V-75) and (IV-79) in recursive form: 

R s ( k |k- ] ) = F ^ T [{k*2)R(k-l|k-2) + v(k-l ) v T ( k-1 ) - HP(k-l j k-2)H T ] 

{ I V-80 ) 


rq s (k|k)r T = ^[(k-l)rQ(k-l |k-l)r T + K(k)v(k)v T (k)K T (k) 
+ P(k|k) - $p(k-l|k-1)* T ] 


(IV-81) 


Sage and Husa devise a suboptimal state estimation algorithm by 
merely using the Kalman filter equations (III-10)-(III-14) and substitut- 
ing the estimates of (1V-80) and (IV-81) for the true values R and rQr T . 
We modify this algorithm in Chapter IV to work for time-varying Q and R. 
The algorithm is given below: 


x ( k | k - 1) = $x(k - 1 Jk - 1) (IV-82) 

P(k}k - 1) = <pP(k - 1 1 k - 1)<? T + rQ(k - Ijk - 1 )r T (IV-83) 

v(k) = y(k) - Hx ( k | k - 1) (IV-84) 

K(k) = P(k | k - 1 ) H T [HP ( k J k - 1)H T + R(k]k - l)]" 1 (IV-85) 

x(kfk) = x ( k ( k - 1) + K(k)v{k) - (IV-86) 

P (' k | k ) = P(kjk - 1) - K(k)HP(k | k - 1) (IV-87) 


R(k + ljk) = pj-yCfk - 2 ) R ( k J k - 1) + v(k)v T (k) - HP(kjk - 1 )H T ] 


(IV-88) 
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rQ( k j k) r T = J-[(k - l)rQ(k - l|k - l)r T + K(k)v(k)v T (k)K T (k) 
+ P { k j k) - ®P(k 1 1 k - 1 )<? T ] 


(IV-89) 


The Method of Magi 11 

The method of Magi 1 1 assumes that the unknown covariances can take 
on a finite number of possible combinations. A bank of parallel station- 
ary Kalman filters is run 3 where each filter assumes a different allowa- 
ble combination of Q and R. The adaptive Kalman filter estimate x(k|k) 
then becomes a weighted sum of the estimates produced by the parallel 
filters [13]. 

We' first assume the system model of { IV-1 )-(IV-4) , where the noise 
covariances are constant and thus denoted as Q and R. The unknown ele- 
ments of Q and R are contained in the vector a; we sometimes use the 
notation Q{ct) and R(a) to indicate that a knowledge of a specifies Q and 
R. Magi 11 assumes that a can take on one of L possible values: 


cte{a-] , c^, • ■ • a^} { IV- 90 ) 

where the ith value has an a priori probability density p[a^]. 

We recall from { I I 1-9) that the optimal minimum variance estimate 
x ( k j k ) is the conditional mean of x(k) given the measurements Y^: 

■ x(k|kj = E{x(k) | Y k ) (IV-91.A) 


= x / xp[x|Y k ]dx 


( IV— 91 ) 


where X is the space of all x(k). Defining A as the space of all a, we 



have: 
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PCX I Y k ] = A f p[x, a|Y k ]da (IV-92) 

By the definition of conditional probability: 

p[x, a|Y k ] = p[x|cc, Y k ]p[a[Y k ] (IV-93) 

Substituting (IV-92) and (IV-93) into (IV-91): 

x(k I k) = x / x A / p[x|a, Y k ]p[a]Y k ]dadx (IV-94.A) 

t 

~ pS ( x / xp[x|a> Y k ]dx}p[ajY k ]da (IV-94) 

where the last step is accomplished by reversing the order of integra- 
tion. From (IV-91) we recognize the term in brackets m (IV-94) as the 
optimal estimate of x(k) given a (Magill calls this the optimal condi- 
tional estimate). (IV-94) thus becomes: 

x(kjk) = A / x(k|k, a)p[a|Y k ]d a - ( IV-95) 

L „ 

= 3 x(k|k, o-OptoJYj • ( I V— 96 ) 

1=1 i i k 

where (IV-96) follows from the quantization of a. 

Magill notes here that the optimal estimation of x(k)' has been fac- 
tored into the linear calculation of a set of conditional estimates and 
the nonlinear calculation of a set of weighting coefficients [13, p. 

434]. The first half of this factorization is easily obtained, since the 
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optimal conditional estimate x(k|k, ct. ) is the linear Kalman filter esti- 


mate produced by assuming ct to be a.: 

x(k|k - 1, a.) = §x(k - Ijk - 1, .a. ) ( I V— 97 ) 

P(k jk - 1, a.) = sP(k - Ijk - 1, a.) + rQ(a.)r T * (IV-98) 

v(k| ai -) = y(k) - Hx( k | k - 1, a.j ) ( IV— 99 ) 

•W(k| ai ) = HP ( k 1 k - 1, ai )H T + R(ct.) (IV-100) 

K(k| ai ) = P(kjk - 1, a.)HV I (k!o i ) (IV-101) 

x(kjk, a-) = x(k|k - 1 , a i ) + K(k|a i )v(k|a i ) (IV-102) 

P(k[k, *.) = P(kjk - 1, a.) - K( k j a- )HP( k| k - 1, a-) (IV-1 03) 


The remaining problem is to find the weighting coefficient p[a. 

We recognize the conditional density p Ca^ | V as the a posteriori 
density of the unknown noise covariance elements given the measurements. 
From Bayes Law we have: 


P [a il Y k ] 


p[Y k ta i ]p[a i ] 

etV 


(IV-104.A) 


= Cp[Y k jct i 3p[cx i ] 


(IV-1 04) 
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where p[Y^] is independent of i and thereby represented by 
whose value is chosen so that the sum of p[a.jY.] over all 

1 K 

now write: 

p[Y k ja i ] = p[y(k), Y r _ 1 Ja.] 

= p[y(k)|Y k _ r a i ]p[Y k _ 1 I a-'] 

We know from ( I V— 1 9 ) that: 

• p[y(k)|Y k _ ^ = N [Hx ( k J k - 1 ) ,' W ( k ) ] 

where x(-k|k - 1) and W(k) are given by the optimal filter, 
a - a- , we have the optimal filter and write: 


■ - p[y(k) |Y k •_ r a.] = N [Hx ( k 1 k - 1,- a.j ) , W(kja-) 

Reapplying ( I V-l 05) to obtain PCY k | cx^ ] for use in (IV-104), 


p[a f | v k ] 


k 

Cp[o-] n p[y(j) ]Y . ,, a,] 

1 j=l 3 " 1 1 


k 

= Cp[a.j n 
1 3=1 


N[Hx( j j j - 1, a t ), WUlcCj)] 


= cp[«.] n | W ( j j a . ) | ^ex p { ryv T ( j | a . ) W” ^ ( j | a . ) v ( j | a • 

1 j=l 1 1 1 1 1 


a constant C 
i is unity. We 

(IV-105.A) 

(IV-105) 

(IV-106) 
Given that 

] ( I V— 1 07 ) 

we have: 

(IV-108.A) 

(IV-108.B) 

)} - (IV-108) 


The adaptive Kalman filter algorithm of Magi 1 1 is now defined. For 
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the ith stationary parallel filter we compute x(kjk, a . ) from ( IV-97 ) - ( I V- 
103) and p[ct.. |Y^] from (IV-108). The adaptive filter estimate x ( k [ k ) fol- 
lows from (IV- 96). ' 

Two modifications will improve the practical implementation of 
Ma gill 1 s algorithm. First, Sims and Lainiotis note that (IV-108) can be 
realized by. a faster recursive form requiring less storage [14]. Vie 
reproduce their result here by rewriting (IV-108): 

p[ ai |Y k ] = Cp[a^] |W(kjct. ) | 2 exp{-iv T (k [cu )W ^ (k Ja^ )v(k|a. ) } 

, (IV-109.A) 

k-1 

X n | W ( j | a - ) ] 2 exp{-lv T ( j [a- )W ^ ( j ja- )v( j |a. ) } 

' * 

_x 

= C [ W ( k ! d - ) I ^exp{-^v^(k )W ^ ( k ja^ ) v( k | a- ) }p[a- | _ -j] (IV-109) 

The second modification results from noting that since the parallel Kalman 
filters are stationary, their gains and covariances reach constant, 
steady-state values. Before actual implementation the gain and covariance 
equations for each parallel filter can be run until K(k[ct..) and W(k|a 1 -) 
reach steady-state values K(a^) and W(a. ). Then ( IV-97)- ( I V-l 03) for the 
adaptive filter can be replaced by: 

x ( k J k - 1, a.) = *x(k - l|k - 1. <» i ) (IV-110) 

v( k ) a^- ) = y(k) - Hx(kj k - 1, a t ) 


(IV-111) 



56 


x(k|k, ot^) = x(k|k - 1, ai ) + K(a 1 )v(k|a 1 ) (IV-112) 

W(a. ) replaces W(k|a^) in {IV-1 09). • ‘ 

The Method of Alspach 

The method of Alspach assumes that the unknown optimal gain K(k) of 
the Kalman filter is a random variable with a posteriori density p[KjY.]. 
Alspach runs a bank of parallel stationary Kalman filters with enough gains 
K- to coverthe space of allowable K. The innovations sequence of the ith 
filter is used, to obtain the density p[K. lY^]. The resulting discretized 
a posteriori density is then used to compute an estimate K(k) of the 
optimal gain for use in an adaptive Kalman filter. 

We assume the system model of (IV-1 ) - ( I V— 4 ) , where the noise covar- 
iances Q and R are unknown constants. For known Q and R the optimal 
state estimate is given by the stationary Kalman filter: 

x(k | k - 1) = $x(k - 1 |k - 1} • (IV-113) 

x(k|k) = x( k | k - 1) + K QpT [y(k) - Hx( k j k - 1)] , (IV-1 14} 

where K^p-j- is the gain K(k) when the Kalman filter equations (III-10)- 

(III-14) are run to steady-state. For the problem here K^p^. is unknown 

and is assumed at time k to be a random variable with an a posteriori 

density p [K [ Y. ] . We next implement a bank of L stationary Kalman 
Kqpt k 

filters running in parallel, where the ith filter has a fixed suboptimal 
gain K.. : 
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:{k| k - 1 , K. ) = $x(k - 1 |k - 1 , K.) 


(IV-115) 


vCkj^) = y(k) - Hx(k|k - 1 , K.) 


(IV-116) 


x(k[ k, K.) = x(k|k - 1, K. ) + K.v(k|K.) - i- : (IV-117) 


The gains (K-j 3 K 2 » • • * K^} are chosen to cover the space of allowable 

Kqp^. We desire to use the observed statistical properties of the ith 

filter to compute the conditional density p., [K. | Y. ] . By computing 

R 0 PT 1 K 

Pi/ [K. |Y,] as a function of K. for a sufficient number of gains we can 
k OPT 1 k 1 

identify the a posteriori density of well enough to estimate its 
value. 

,He define W^p-j. as the steady-state innovations covariance of the 
optimal filter. Alspach first solves for the joint a posteriori density 
of Kqp^. and Wqpjj which by Bayes Law becomes: 


Pi/ 

n opt 


’ W 0PT 


[K, W 



p[Y. 1 K 3 W]p K w 
K K 0pT , W QpT 

p[VJ 


[K, W] 


(IV-118.A) 


.= Cp[Y k |K, W]p[K, W] (IV-11S) 

where p [ Yj^] is constant for all K, W and therefore replaced by C. 

Pv w [K, W] is the a priori density of K and W, representing any 
^OPT’ W 0PT 

knowledge of and Wgp-j. without measurement information. The sub- 
scripts on p[*] are dropped where no confusion results. From the defini- 
tion of conditional probability we can rewrite {IV-113): 
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p[K> W|Y k ] = Cp[Y k |K, W]p[W|K]p[K] . (IV-119) 

We then write: 

PCYjjK, W] = p[y(k), Y r _ ^ | K, W] (IV-120.A) 

a pCy(k) |Y k _ r K, W]p[Y k _ 1 |KW] (IV-120) 

From ( I V— 1 9) we know that: 


P[y(k)jY k _ ^ = N[Hx( k J k - 1, K 0pT ), W QpT ] (IV-121) 

Given thax Kg p -p = K, = W, we can therefore impute: 

p[y(k)|Y k _ r K, W] = N[Hx(kj k - 1, K) , W] (IV-122) 


Applying this result to (IV-120) we have: 


pTy.Ik, w] =. n p[y(j)|Y, , , k, 

j=l J 


W] 


(IV-123.A) 


= n N[Hx{ j | j -1, K)J] 
j=l 


( IV- 1 23 ) 


We now rewrite (IV-119) 


PtX W,|Y k ] = Cp[W | K]p[K] _n |W| 2 exp{-|[y(j) - Hx(j[j - 1, k)] T 


j=l 


(IV-124.A) 


,-l 


x W [y( j) - Hx( j j j - 1, K)]} 
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k 

= Cp[W|K]p[K]|W| 2 exp{4 2 v T (j|K)ir 1 v(j|K}} (IV-124) 

Alspach now specializes the adaptive Kalman filtering development 
to the case of scalar noise and measurements. This is general enough to 
cover our specific aircraft model of (III-l )-( II 1-4) . Q, R, and Wq^ are 
now scalars, and (IV-124) becomes: 


* k k 2 

p[K, W|Y k ] = pCW)K]p[K]W'2exp{-l ^ 


(IV-125) 


We define the sample covariance of v(kjK) by: 


^ 1 * o 

W(kjK) - ~ I v^{ j i K) 
k j=1 


{ IV- 1 26 ) 


For stationary Gaussian noise W(kjK) is an unbiased estimate of W. (IV- 
125) thus becomes: 


k 

p(X WjY k ] * p[W|K]p[K]WTexp{-|^i^-} 


(IV-127) 


We now express W opT in terms of K QpT by using (III-32): 


“opt = R/[1 - HI W 


(IV-128) 


Assume that the only a priori information we have about Q and R is that 
they are bounded by the values and Rj^. T^ en given the only 

information we have about hgpy is that it lies somewhere between zero and 
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W MAX^ K 0PT^ " WP “ HK 0PT^ 


( I V— 129) 


Alspach therefore assumes that the conditional density of Wg pT given K 0 p T 
is uniform: 


p(W|K) = 


ff : 0 _< “ ; w max (k) 

MAX mA 

0: otherwise 


(IV-130) 


Substituting (IV-130) for p[K|W], we can obtain p[K|Y^J by integrating 
the joint conditional density of (IV-127) over the range of W: 


p[K[Y k ].= / p[K, W|Y k ]dW 


(IV-131 .A) 


W r MAX^ 


Cp[Kj 

w max ( k) 


W( k I K) . 


W 2exp{ }dW 


(IV-131) 


Us 


1/ A 

ing the variable of Integration z = ^(kjKj/W, (IV-131) becomes: 


p[K Y k ] = ^l)-[W(klK)r^ k ‘ 2)72 j z" (k " 4)/2 e" Z dz (IV- 132. 
'•I AX 


A) 


= ^ [ | k w(k|K)]- (k " 2)72 WT(k s K) 


MAX 


(IV-132) 
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where 

. a = ■|w(k|K)/W f|AX (K) ( IV- 1 33) 

WT(k, K) is the integral of (IV-132.A) and is a function of the ratio 
W(kjK)/W MA ^(K). We define M as the value (k - 4)/2. When H is an inte- 
ger (meaning k is even and greater than 2) we can evaluate WT : 

WT(k, K) = / z M e Z dz (IV-134.A) 

a 

M ,H - j 

= M!e" a I jyj (IV-134.B) 

j=0 v " JJ ’ 



- )!ex P { -|w 


W(k 

MAX 



(k-4)/2 

z 

j=0 


r k j(k, K) , (t 

L 2H^IkT J 

l }L f A - J]! 


41/2 - J 


(IV-134) 


Alspach plots WT(kjK) as a function of W(k|K)/W^^(K) for different 
values of k. His results are reproduced in Figure IV-2. He notes that 
for k above lOOOjWTCk, K) can be approximated as a unit step which falls 
to zero at W(k|K) = W^^(K) (This is not a bad assumption even for 
k = 50 or 100). Thus WT(k, K) acts to discriminate against gains for 
which the sample innovations covariance exceeds the maximum value: 


p[K|Y k ] = | 


i®b tS ' k iwr (k ' 2)/2 = hihw .<»„(« 


0: W(k|K) > W mx (K) 


( I V— 1 35) 
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We now assume no a priori information about Q and R 


and become infinite. Since nothing is known about K 

made uniform over the limits of allov/able gain. W^^(K} 

A 

infinity, so that WT(kjK) becomes unity for all W(kjK). 


, so that Q mx 
0 pT s p[K] can be 
now approaches 
We now have: 


p[K|Y k ] = 


|c[W(k|K).] ^ allowable gain 


0, otherwise 


(IV-136) 


The adaptive Kalman filter of Alspach is now ready for implementa- 
tion. For each of the parallel suboptimal filters the estimate 

A 

x(k]k, K^) and innovations residual v( k J K - ) are obtained from (IV-115)- 
(IV-1T7). The a posteriori gain density p[K.. {Y^] is then computed by 
(IV-136), using the sample innovations covariance of (IV-126), which 
Alspach writes in recursive form: 


W(k|K) = l[{k - 1 )W(k - 1 JK) + v 2 (k|K)] (IV-137) 


The adaptive filter implementation’ is greatly simplified when we 
recall from Chapter III that for scalar noise and measurements the opti- 
mal gain Kgp-j. is known when only. one element is specified, (^(k) is 
given as a function of K-j(k) in (II 1-42 ) for our aircraft model). Also, 
the gain Kgp-j- is always bounded. We can therefore implement the parallel 
filters by uniformly .incrementing the first gain element between limits 
K-. and K-, , allowing K..('2 < j < n) to be determined by the appro- 

priate functional relation. We can therefore use (IV-136) to find 
p[K, | Y. ] and thus estimate K, . Alspach does not specify what 

i K 'opt 
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estimation scheme to use in computing K^(k). However, since we know the 
conditional mean £{K] |Y^,} to be the optimal minimum variance estimate 
from (I 1-9), we use it here: 


K-j(k) 


L 

t K.pIXfY.] 
i=l 1 1 K 


(IV-138) 


A A 

K-j(k) automatically defines K(k), which is then used in the adaptive 
filter of (IV-5) and (IV-6). - 

Alspach admits that his algorithm may be impractical for use in a 
general-purpose digital computer where the L suboptimal stationary fil- 
ters must be implemented serially. However, in a special-purpose paral- 
lel implementation the stationary filters can run simultaneously, produc- 
ing a fast adaptive algorithm. He also notes that, though similar to the 
parallel filters method of Magi'll, his algorithm is simpler, requiring 
fewer parallel paths. Consider the scalar noise case where Q and R are 
both unknown,. For n possible values of Q and m values of R, we would 
need n x m parallel filters in Magi 1 1 1 s algorithm. The number of allowa- 


ble Q and R values may increase further if we do not know their upper 
bounds and Rj^, which can approach infinity. In Alspach 's algo- 
rithm the only unknown is K^ , which is always bounded. -We need only to 
use enough parallel filters to adequately cover the range of allowable K-j 
values. 


The Minimum Innovations Covariance Method 

The last adaptive Kalman filtering method presented here- is an’ 
intuitive scheme which has some theoretical backing. As in the methods 
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of Magi 11 and Alspach, a bank of parallel fixed-gain Kalman filters is 
run. The gain of that filter with the minimum innovations sample covari- 
ance is chosen as the gain K(k) for use in the adaptive Kalman filter. 

We again assume the system model of (IV-1) to (IV-4), but with 
scalar noise and measurements, so that Q and R are unknown scalar con- 
stants. From (II 1-13) we recall that the steady-state innovations covar- 
iance W is minimum when the suboptimal filter gain becomes the optimal 
Kalman gain. Let us revisit Figure III-l and IV-1, where P-j.|(k|k); the 

(V 

mean-square error in e(k|k), and W are plotted as functions of suboptimal 
>\ ' 

gain K-j for the aircraft system of (III-l )-(III-4) with stationary noise 
(K, is given by (II I -42 ) ) . Not only are P-j-j(kjk) and W minimum for 

A A 

K-, = K, , but both increase monotonically when K. either increases or 

1 ‘opt 1 

decreases from the optimal gain. We can assert that the innovations 
covariance W is a direct indicator of a suboptimal filter's error 
performance. 

We now implement a bank of parallel fixed-gain Kalman filters, 

where, as in the algorithm of Alspach, the gain K-j is incremented between 

the limits K, and K-, . The ith filter is realized, just as in 

MIN ‘MAX 

Alspach, by (IV-115)-(IV-11 7) , The sample covariance W(kja^ ) of the 
innovations sequence is computed from (IV-137). We note that W(k|K) is 
the standard covariance estimate for a stationary scalar process with 
zero mean, and is therefore our best estimate of U(k[K). It would thus 
seem reasonable to choose that parallel filter with the lowest value of 
W(k|K. ) as the one whose gain is closest to K, . We therefore choose 

1 "opt 

the gain of this filter as the gain K(k) of our adaptive Kalman filter: 
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K(k) = K-: W(k|K t ) < W ( k j K d ) , 1 < j < L (IV-139) 

The adaptive filter is given by ( I V-5 ) and (IV-6). 

This intuitive scheme has a theoretical appeal when we recall (IV- 
136), which gives p[K^ | for the ith parallel fixed-gain filter. The 
gain K. which minimizes W(k|K. ) is the same gain which maximizes the a 
posteriori density pfK^. [Y^] } as derived by Alspach. Instead of the con- 
ditional mean estimate of Kq P -j., we are choosing the maximum a posteriori 
estimate of Kgpj- The intuitive, minimum sample covariance estimate K(k) 
may therefore be considered the MAP gain estimate of Alspach. Of course, 
this MAP- estimate does not require the calculation of pfK^ |Y.j, and is 
therefore easier to implement than the method of Alspach. 



CHAPTER V 


ADAPTIVE KALMAN FILTERING WITH TIME-VARYING NOISE STATISTICS 


In Chapter IV we presented four methods of adaptive Kalman filter- 
ing for use in problems where the noise is stationary. We recall our 
general system model of ( I V- 1 )- ( IV-4) . The noise covariances Q(k - 1) 
and R(k) were unknown' constants. In this chapter we remove the station- 
ary noise assumption and allow Q(k - 1) and R(k) to vary with time. We 
modify the adaptive filtering methods to work for the non-stationary 
noise case and then specialize them to our aircraft system model. 

. The Method of Alspach ' 

We recall that Alspach assumes that the measurements and noise 
terms of the system model (IV-1 )- { IV-4) are scalars. The a posteriori 
density of the optimal gain K(k) is found by computing the sample innova- 

A 

tions covariance W(k|K. ) for each of the parallel stationary filters: 

W(k|K.) =^[(k - l)W(k - 1|K.) + v 2 {k|K.)3 (V-l) 

Alspach points out that as k becomes large, the present innovations 
residual v(k|K^) has little effect upon the value W(k|K^). In order to 
prevent W{k[K.) from becoming insensitive to new information he suggests 
the following change [15, p. 553]: 
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W(k[Kj) 


l[(k - l)W(k - 1 | K . ) + v 2 (k|K.)]; k < N 
1[(N - l)W(k - 1|K.) + v 2 (k[K-)]; k > N 


(V-2) 


where N is chosen for the stationary noise case such that W(kjK.) is ■ 
within some acceptable r.m.s. deviation of the true covariance W(k|K. ) 
for k > N. We can view (V-2) as a fading-memory estimate of W(k|K. ), 
where old innovations residuals are deweighted: values v(j|K. ) will have 

little effect on W ( k J K^. ) for j < (k - N). W(k|K..) becomes the output of 

a first-order lowpass filter with input v (k) and time constant NAt. 

If the noise covariances Q(k - 1) and R{k) are slowly changing with 
time, we can approximate them as being constant for N iterations. We can 
then estimate the state with stationary noise methods where only the last 
N innovations residuals are used. Alspach has done this by using (V-2) 
to estimate W(k|K. ) as it changes with Q(k - 1) and R(k). Of course, the 
more slowly changing the noise covariances are, the larger N becomes, ' 
making W(kjK.) more accurate.' Alspach also modifies (IV-136) for comput- 
ing the a posteriori density of the gain: 




C[W{k|K.)]" (k ’ 2)/2 , k < N 
C[W(k|K. )] -(N “ 2)72 , k > N 


(V-3) 


The adaptive algorithm of Alspach remains the same as for the sta- 
tionary case, with (V-2) replacing (IV-137) and (V-3) replacing (IV-136). 
Alspach adds another modification by restricting the range of W(k|Kj ) 
among parallel filters. This is done to enhance the adaptive filter's 
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ability to follow ti me changes in the noise covariances, and is illus- 
trated by example: 

Consider two filters in the parallel filter bank, one with very low 
gain K. and the second with high gain K, . Consider also the case where 

V# ' 

the ratio Q(k - 1)/R(k) is large, so that the optimal gain sis near K, . 
Recalling Figures III-l and IV-1 , we expect W(kjK^) to be low while 
W(kjK.) becomes very high, indicating a diverging filter. p[K.}Y.] will 
be nearly zero, so that only the higher gains contribute to'K(k). Now 
assume that Q(k - 1)/R(k) suddenly becomes small, so that the optimal 
gain is near K.. The true innovations covariance W(k[K.) becomes small, 

J J 

but the sample covariance W(kjK.) will not show this effect for a consid- 

J 

erable time; the old residuals v(k|K.) taken while the filter was diver- 

vJ 

gent must be deweighted and replaced by new residuals of lower covari- 
ance. Such a process could require more than a time constant of the 
fading memory filter of (V-2). 

A A 

Alspach has therefore placed a ceiling on W(kjK-), If W(k|K, ) is 

J * 

A 

the minimum sample covariance among all parallel filters, then W(kjK-) is 

J 

not aliovied to exceed an upper bound f^^ x W(k)K-j), Whenever this limit 

A A 

is exceeded, we replace the estimate x ( k j k , K.) with x(k|k, K, ) . This 

d 

modification allows p[K-|Y. J to quickly become significant when the opti- 

J K 

mal gain suddenly shifts toward Kj. 

Alspach also modifies his algorithm to allow the value N to adapt 
to changes in the time variations of Q(k - 1) and R(k). A fading memory 
estimate W(kjK) of the adaptive filter's innovations covariance is com- 
puted using (V-2). A second estimate W^kjK) is computed by replacing N 
.in (V-2). with a smaller time constant Ng (we could make some fraction. 
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say 20 percent, of the nominal value of N). If N is large and 
Q ( k - 1)/R(k) suddenly changes, the small time constant filter will soon 
detect this change by changing ( k J K) . When W(kjK) and ( k } K) differ 
by more than an allowable amount, N will be decreased. 

Alspach uses the following procedure for changing N: for station- 

ary noise we know that the variance in the unbiased estimate W 2 (k]K) is: 


' Var{W 2 (k|K)} * 2W 2 (kjK)/N 2 . 

(V-4) 

Assuming W(k|K)*to be our best estimate of W(k|K), we 

can measure stand- 

ard deviation in '^(kjK): 


a 2 = W(k|K)/r1y2 

( V— 5) 

Alspach then modifies N according to the rule: 


A = J W ( k | K ) - W 2 (lc|K) | 

(V-6.A) 

- IF A < ^2' N -»• N + N 2 

(V-6.B) 

IF A > 2a 7 : N -* N - Integer[^-Nj 

(V-6.C) 


The above procedure works for situations where changes in Q(k - 1) and 
R(k) produce a wide dynamic range in the values of W(k|K). In the simu- 
lation testing described in Chapter VI this was not the case, .with W(k|,K) 
never changing by more than 25 percent. N was chosen experimentally from 
various simulation runs and left constant. 



Figure V-1 is a block diagram of the adaptive algorithm of Alspach 
for the specific MLS aircraft model o.f (III-l ) - ( 1 1 1-4) . The algorithm is 
implemented as a computer subroutine, where during the kth iteration the 

* m 

measurement y is received and estimates 0 ^ and.0^ are computed and 
returned to the main program (the subroutine is not given the- noise 
covariances Q(k - 1) and R(k)). The adaptive filter has L parallel sta- 
tionary filters, implemented according to (IV-115)-(IV-117) . The condi- 
tional density p[K-[Y,,] is computed for each filter according to (V-3), 
using the sample innovations covariance W(k|K.) of (V-2). The adaptive 

' A 

Kalman filter is then updated using the gain K computed from (IV-138). 

The block diagram shown here is for an algorithm using constant N. 
FIRST is a logical variable which is TRUE until k is greater than N. We 
should point out that since this algorithm will be implemented as a sub- 
routine in a digital computer simulation, the parallel filters. must be 
run serially. In an actual parallel implementation, the-i— loops in the 
block diagram would not exist: the parallel filters would run 

simultaneously. 

We now procede to modify the other three adaptive filtering algo- 
rithms for use in problems where the noise statistics are time varying. 

We make use of Alspach's fading memory approach to deweighting old inno- 
vations information. Each algorithm is then specialized to our aircraft 
problem model of (III-l )-( III-4) . 

The Minimum Innovations Covariance Method 

As stated in Chapter IV, this method is actually that of Alspach. 
where the maximum a posteriori (MAP) estimate of the gain K(k) is used by 




72 



Figure V 


Adaptive 'Algorithm of Alspach 








Figure V - l.A Algorithm of Alspach (continued) 
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the adaptive filter. The adaptive filter algorithm looks like that of 
Alspach, except that the conditional density p[K.. JY^] is not computed. 
Instead, the adaptive gain K is just set equal to the gain' of the paral- 
lel filter with lowest sample covariance W(k|K.). We can therefore use 
the block diagram of Alspach in Figure V-l , where the only modifications 
necessary are between points A and B. These modifications are shown in 
Figure V-2. 

The Method of Sage and Husa 

The suboptimal algorithm of Sage and Husa is that of a discrete 
Kalman filter where the unknown noise covariances Q(k - 1) and R(k) are 
replaced by the estimates of (IV-75) and (IV-79). We can make these 
estimates responsive to changes in the noise covariances by using only 
the last N innovations residuals: 


rQ s (kjk)r T 


~ Z .K(j)v(j)v T (j)K T (j) + P{j|j) - »P(j-l |j-l)® T 
" j=k-N 


(V-7) 


k-1 

R,( k|k - 1) = l z v(j)v T (j) - HP ( j I j - 1)H T 
5 N j=k-N-l 


( V— 8 ) 


We now use recursive approximations : 


rQ s ( k J k ) r T 


£[(k - 1 ).rQ s ( k - l[k - l)r T + K(k)v(k} v T (k)K T (k) 

+ P(k)k) - <&P ( k - 1 J k - 1)* T 3. k < N 

jj-[(N - 1 )rQ s (k| k)r T + K(k)v(k)v T (k)K T (k) 

+ P(k | k) - aP(k - 1 |k - l)i> T L k > N 


(V-9) 




Figure V - 2. Modifications to Block Diagram of Alspach for 

Inplementation of Minimum Innovations Covariance 
Filter 
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]~ T [(k-2)R s (k-l | k— 2 ) +v ( k ) v T (k)-HP(kj k-l )H 1 ] » k < N 
^-[(N-l )R S ( k-1 |k-2)+v(k)v T (k)-HP(k|k-l)H T ], k > N 


R_{k|k-l) = < 


(V-l o) 


We note that the equations in (V-9) and ( V-l 0 ) for k < N are the recur- 
sive equations (IV-88) and (IV-89) in the original algorithm of Sage and 
Husa. We can obtain a more practical form of (V-9), in terms of quanti- 
ties already computed by the Kalman filter, by recalling from (IV-83) 
that $P(k - 1 1 k - l)® 1 * equals P(kjk - 1) - rQ s ( k - 1 J k - l)r T : 

y 

rQ s (k - Ijk - l)r T + ±-[K(k)v(k)v T (k)K T (k) 

T + P(kjk) - P(kjk - 1)], k < ,N 

rQ s (kjk)r ~ . (V-Tl ) 

rQ s (k - l|k - 1 )r T + l[K(k)v(k)v T (.k)K T (k) 

+ P ( k | k ) - P(k|k - 1)], k > N 

The adaptive algorithm of Sage and Husa is now given by the origi- 
nal algorithm of ( I V-82) - (IV-89) with (V-10) replacing (IV-88) and (V-ll) 
replacing (IV-89). Figure V-3 is a block diagram of this algorithm for 
the specific aircraft model of (III-l ) - ( I I I -4) . The Kalman filter equa- 
tions (IV-83)-(IV-87) are given by the specific equations (III-18)— (III- 
28), with- rQ s (kjk)r replacing rQ(k)r . R( k) is known in this problem, 
and so R s (k|k) is not needed. 


The Method of Magi 1 1 

Here we modify the algorithm of Magi 11 for the case of scalar mea- 
surements and noise. We make the conditional density p[a- | Y^3 in 
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P00 ^ quality 



Figure V - 3 . Adaptive Algorithm of Sage - Husa 
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(IV-108) a function of only the N most recent innovations residuals: 


p[a,-|y 


- Cp[cj-j] 



(V-12) 


We then use the recursive approximation: 


phiV - 


1 i. v (k|a. } ' 

CW-2(k|c,,)e X p {T - ifrq J 7 }p[a i |Y k . ,] 

, 2,., , M - 1 

±]_ -I v (k|a*) m 

CW 2 ( k | ai 0 e xp4 w ^l ] -}(p[a.|Y k _ ,3) 


, k < N 


, k > N 


(V-1 3) 


where p[a.jY 0 ] = p[o f ]. 

We recall that for our specific problem model (III-l )— (III-4) only 
the plane -noise covariance Q(k - 1) is unknown. We can therefore set the 
unknown parameter vector a equal to the scalar Q(k - 1). We implement L 
parallel Kalman filters: the ith filter uses the true value R(k) and an 

estimate for Q(k - 1). 

The ith parallel filter is not stationary, since R(k) is time vary- 
ing. The gain K - ( k ) is not a steady-state value as in Chapter IV, and 
this requires running all the Kalman gain and covariance equations. This 
problem can be avoided in our case by noting that is never explicitly 
used in computing p[Q. |Y,J; only K(k)Q.) and W(kjQ.) are needed (here we 
have replaced ct. with Q. ). Since for every Q. and R(k) there exists a 
unique K(kjQ.), we might ask why the gain could not be a conditioning’ 
variable instead of Q. . We embrace this approach here. For our problem 
K^(k) is bounded between 0 and 1. We therefore implement a bank of 
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parallel stationary filters with first gain elements K, uniformly spaced 

i 

between these limits. K 2 is given by (1 11-42) . We can use ( II 1-32) to 

•i 

•obtain W(k [ K. ) : 

W(k|K.) = R(k)/[1 - ^ ] { V-l 4) 

The conditional estimates for the ith filter are given by: 
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equations must be run for the ith parallel filter, since the gain K^- is 
fixed. Also, as noted by Alspach, the optimal gain is -bounded, while 
Q{k - 1) may take on any positive value. For the general problem where 
Q(k - 1} and R(k) are both unknown, this modification of Magill cannot be 
used. Given only K. , we do not know W(k|K.) and thus p[K*jY^] remains 
unknown. We must then resort to using parallel filters where various 
combinations of Q and R are assumed. 

Figure V-4 provides a flowchart of the modified Magill "algorithm 
for the aircraft system model (III-l )-(III-4) . As for the other algo- 
rithms, the adaptive filter is implemented as a subroutine, receiving 
y(k) and R(k) and returning estimates Q^p(k k) and g( k [ k ) . The a priori 
density of the gain K-j(k) is assumed uniform. 








Figure V - 4. Adaptive Algorithm of Magi 11 






Figure V - 4. A Algorithm of Magill (continued) 




CHAPTER VI 


COMPUTER SIMULATION TESTING 


We now develop and employ a digital computer simulation' for testing 
the candidate adaptive filters for our aircraft landing problem. We rea- 
lize at the onset that poor performance can occur for one of two reasons: 
first, an adaptive Kalman filter may be a poor estimator, given the sto- 
chastic state model for which it was developed; secondly, the assumed 
state model may inadequately describe the physical system for which the 
adaptive filter is used. Our testing is therefore conducted in two 
phases^ We first simulate the state model (III-l ) - ( 1 1 1 -4) , repeated 
bel ow : 


| e( k ) j p At 

Le<k)J i_P 1, 


e(k - 1) 


0 




_0(k - 1)_ 




m(k - 1) . (VI-1) 


y C k } = [ i o ] 


~e(k)~ 

'e(k) 


+ v(k) 


( VI — 2 ) 


• -P[o)(k - 1)] = WN[0, Q(k - 1)] (VI-3) 

p[v(k}] = WN[0, 'R(k)] (VI-4) 

At each new time increment the measurement y(k) and the error covariance 
R(k) are sent to the candidate adaptive filter, which computes state 
estimates ©(kjk) and o(kjk). In the second test phase we remove the 
state model (VI-1) and update e(k) deterministically, as in (II-l): 

83 
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0 s f (k) (VI-5) 

where f(*) describes the ‘evolution of 0(k) for aircraft motion along a 
given flightpath. The measurement model (VI-2), (VI-4) is retained, and 
the adaptive filter is retested. In this two-prong approach we establish 
the performance of each candidate filter for both' the assumed state model 
and the actual landing approach. 

In simulating the stochastic system model (VI-1 )-(VI-4) we must 
select realistic functions for Q(k) and R(k). R(k) is the covariance of 
the error in the estimate- y(k) computed by the locally optimum estimation 
algorithm in the envelope processor. This covariance has been computed 
as a function of receiver signal -to-noise ratio in earlier simulations 
[5, pp. 25-2?]. Q(k) is the covariance of the Gaussian white noise driv- 
ing the system. We recall from Chapter II that our state model (VI-1) 

»• 

was derived from a continuous-time model where acceleration e(t) was rep- 
resented. by white noise. The mean- square of the noise was set equal to 
the square of the acceleration. From (11-11) and (11-22) the discrete 
model noise covariance is given by: 

Q(k - 1) = 4te 2 (t k ) (VI-6) 

Of course, the aircraft does not know ©(t^), so Q(k - 1) is unknown as 
well . 

In light of (VI-6) we use the following scheme for propagating the 
state model of (VI-l). The true acceleration e(t k ) is computed in a sub- 
routine for a typical deterministic flightpath and passed to the main 
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program. Here Q(k - 1) is computed using (VI-6). The system state is 
then updated using, a sample from, a white Gaussian population with covari- 
ance Q(k - 1). We are careful here not to confuse ©( k) and q( k ) with 
©(.t^.). The former are states of a stochastic process driven \hy white 
noise, ©(t^), a deterministic quantity, is a tool for setting Q ( k - 1) 
in the simulation and has nothing to do with the state. , 

The same subroutine which computes ©(t^) for aircraft motion along 
a given flightpath also computes ©{ t ^ ) and ©(t^J. . In the second phase of 
simulation., where ©(k) is updaied deterministically, these values are 
merely passed to the main program, which sets ©(k) equal to e(t^) and 
e(k) equal' to s(t,). We now address the task of realizing a suitable 
flightoath for updating both Q(k) in the stochastic case and e(k) in the 
determinisxic case. 

The Landing Approach 

Before assuming a test landing pattern, we first place some 
restrictions on the set of allowable flightpaths. Let f(k) describe the 
evolution of e(k) as the aircraft travels a given flightpath. We recall 
from Chapter II that, while unknown, f(k) is a member of a known class of 
functions. We restrict this class to include those functions attributab- 
le to aircraft motion along a restricted family of flightpaths. This 
family includes what we assume to be reasonable flightpaths, thus ruling 
out unrealistic approaches for which adaptive filtering could not work. 
For example, a missed approach where the aircraft crosses the runway at 
high speed within a mile of the azimuth antenna, produces very high and 
rapidly changing values of ©(t). The resulting covariance Q(k) in the 
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system model would be too rapidly varying to be followed by an adaptive 
filter. We therefore place the following restrictions on the landing 
approach: , ' 

1. Maximum airspeed = 200 knots 

2. Minimum turn radius = 1 N. mile 

3. Flightpath must be coincident with runway 
centerline before runway is reached (no 
missed approaches) 

A 

We assume these conditions to be those of a worst-case approach. 

Given these restrictions, e(t) has been observed in simulation to 

2 -r 

remain below 0.1°/sec. [6, p. 40]. From (VI-5) we can therefore place 

an upper limit on Q(k): 

= (.075)(.l) 2 = 7.5 x 10' 4 (VI-7) 

where At = .075 seconds, the time between the start of successive azimuth 
scans (at a 13 1/3 Hz update rate). 

We now return to the task of finding a suitable test flightpath 
within the above restrictions . Figure VI-1 shows a representative land- 
ing approach selected for this simulation. The aircraft travels at 120 
knots along an S-curve flightpath, staying on runway centerline for the 
last 3 N. miles of the approach. The runway is 2 N. miles long, with the 
azimuth antenna at the stop end. 

We have developed a FORTRAN computer subroutine for computing 0(t) 
and its derivatives as the aircraft follows an S-curve approach of varia- 
ble dimension. This general flightpath, shown in Figure VI-2, has the 
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Figure VI - 2. General S - Curve Flightpath. 
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following programmable parameters: 

v .• = airspeed (knots) 

d.j , d 2 , d 3 = lengths of the straight segments of 
the flightpath, as shown in Figure 
VI-2 (N. miles) 

r^, = turn radii (N, miles) 

1 = distance (N. miles) from azimuth 

antenna to that point where the 
approach first coincides with the 
runway centerline (5 N. miles for ■ 

Figure VI-1 ) . 

A flowchart of the flightpath subroutine is given .in Figure VI-3. 
The various straight and curved sections of the approach are labeled from 
A to E on. both the flowchart and Figure VI-2. The simulated flight runs 
from time t^ to t^, with t^ through t^ marking transition times from one 
flightpath section to the next. The subroutine receives the present azi- 
muth scan number and computes the time t. Based upon which flightpath 
section the aircraft is currently following, its cartesian coordinates x 

i 

and y and their derivatives are computed. The subroutine then computes 
6(t) , 0(t) , and @(t) : 


G = arctan(y)x) 

(VI-8) 

0 - (xy - xy)/(x 2 + y 2 ) 

(VI-9) 

- xy - 20 (xx + yy)]/(x 2 + y 2 ) 

(VI-10) 
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The values ©(t), o(t), and e(t) are returned to the main simulation pro- 
gram, along with the logical variable FINAL, which is set to TRUE when 
the simulated approach has been completed. 

Figure VI-4 shows the time functions ©(t), ©(t), and e(t) for the 
S-curve flightpath of Figure VI-I . The simulated approach begins with 
the aircraft 2 N. miles ahead of the first 90° turn and ends 2 N. miles 
beyond the second turn, or 1 N. mile before the runway is reached 
(d^ - d^ = 2 N. miles). At 120 knots the aircraft covers 8.88 N. miles 
in 265 seconds, or 3550 azimuth scan periods (scan update rate = 

13 1/3 Hz). Figure VI-4 also shows the plant noise covariance Q ( k ) for 
the stochastic modeling of this flightpath, computed from (VI-5). 

Computer Simulation Structure 

We now- describe the actual test simulation, implemented as a 
FORTRAN computer program. A flowchart of the overall simulation is given 
in Figure VI-5. At the kth scan period the main program calls the 
flightpath subroutine, which updates ©(t) , e(t), and ©(t) according to 
deterministic aircraft motion along the S-curve flightpath. The main 
program then updates the state values ©(k) and ©(k), either stochastic- 
ally with white noise or deterministically , depending on the value of the 
logical variable MODEL. When MODEL is TRUE, the state is updated with 
noise according to (VI-1). Q(k - 1) is computed from ©(t) using (VI-6). 
The white noise term w(k - 1) is then obtained from Q(k - 1) and the 
output of GAUSS, a subroutine which uses the machine random number gener- 
ator to produce independent samples of a standard Gaussian population 
(zero mean, unity covariance). When MODEL is FALSE, the state is updated 
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deterministi cal ly by merely setting e(k) to e(t) and e(k) to e(t) (x-j and 
*2 are used for ©( k ) and 0 (k) in the actual program to avoid confusion 
with o(t) and o(t) when the stochastic model is used — see flowchart in 
‘Figure VI-4). ' 

The main program now uses GAUSS and R(k) to produce an additively 

f* 

corrupted measurement y{k) of the state* as -in (VI-2). Q(k - 1), R(k), 
and y(k) are then sent to the Kalman filter subroutine., which computes ‘ 
the optimal state estimates 0 Qp^.(k|k) and 0 Q p ^(k|k) (these estimates are 
optimal when the assumed state variable model (VI-1 )-( VI— 4) is correct). 
R(k). and y(k) are then sent to the candidate adaptive Kalman filter sub- 
routine r which, without knowledge of Q(k - 1), computes the suboptimal 

- A A . 

state estimates 3 ^, D (k|k) and e^ D (k|k). The main program computes errors 
in estimates of e(k) and increments k to the next scan period. 


Simulation Testing: The Stochastic Model Case 

The S-curve flightpath of Figure VI-1 is used in both the stochas- 
tic and deterministic phases of simulation testing. Here we use Q(k) . 
from Figure VI-4 to update the state variable model (VI-1). 

As a result of the restrictions placed on the family of allowable 


flightpaths, we can limit the adaptive gain K-j . We recall from ( I V— 7) 

— 4- ' 

that Q(k) is bounded at = 7.5 x 10 . We also recall from Chapter 

III that the Kalman gain K-j(k) is a monotone increasing function of' 

Q(k - 1)/R(k), Given Q^, we can thus limit K-j by placing a lower bound 
on ?*. In earlier simulations we have found the r.m.s. error in the enve- 


lope processor estimate y(k) to remain above .01° for the expected range 
of signal -tc-noise ratios (20 db or less) [5, p. 27]. We-. assume a lower 



limit on r.m.s. error in y(k) of .005°, yielding the bound 

r 

RfilN = 2-5 x 10 * Whefl ^MAX and R MIN are used ln 0Ur s ^ stem mode ' 1 tlie 
Kalman filter has a steady-state gain K-j = .602. We assume the upper 

bound K, = .625. 

*MAX 

Three of the adaptive filters tested here use a bank- of parallel 
stationary Kalman filters. In each case we use 24 parallel filters, 
incrementing K^ uniformly from .05 to .625. For each of the three adap- 
tive filters we use .a fading memory time constant N of 80 scan periods in 
computing the sample innovations covariances. This value' has been chosen 
experimentally by studying the effect of different values of N on filter 
performance for various flightpaths (from typical to worst-case 
approaches). . The adaptive filter of Sage and Husa uses a constant N of 
30 scan periods. ' 

In-addition to the adaptive filters of Chapters VI and V we also 
test a suboptimal filter which does not adapt'to chang es in Q(k). This 
estimator is merely a Kalman filter which uses the true value R(k),but 
which replaces the unknown Q( k) with the limit of 7.5 x 1(T^, from 
(VI-7). Such an estimator, which is much simpler than an adaptive 
filter, has some intuitive appeal. Since R(k) is known, the filter gain 
K-j(k) is always greater than or equal -to the optimal gain K-j ( k ) (as 
Q[MX - Recalling Figure III-l, such a filter, while not'always 

optimal, always has a mean- square error in © ( k j k ) less than R(k). We test 
second constant-Q filter which has knowledge of maximum acceleration e(t) 

for the actual flightpath to be used. For this simulation the S-curve 

•• o 

approach produces a peak o(t) just above .01°/sec , while Q(k) peaks at 
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8 x 10 (see Figure VI-4.C, D). We ha\e assumed in the problem defini- 
tion of Chapter II that nothing is known about the actual flightpath 
except that it belongs to a given family of approaches. We nevertheless 
include this filter for test comparison with the other candidate filters. 

We assume a constant R(k) of 10"^ for the simulation (r.m.s. error 
in y(k) = .01°). This is not too realistic, because the si gnal -to-noise 
ratio slowly rises as the aircraft approaches the runway. A slowly 
decreasing function for R(k) would seem more reasonable. But judging 
from Figure VI-4.D, we would expect the more rapid variations in Q(k) to 
cause the most difficulty in adaptive estimation. A constant-R simula- 
tion should give a fair indication as to whether or not adaptive filter- 
ing will work. 

The candidate adaptive filters have been tested in a FORTRAN simu- 
lation on a PDP-1103 computer. The stochastic state model (VI-1 )- (VI-4) 

is implemented, with Q(k) updated as shown in Figure VI-4.D. R(k) is 
-4 

constant at 10 . The nonadaptive filter using Q ^ has a steady-state 

gain K-j = .476. The second constant-Q filter, which has knowledge of 
maximum acceleration for the actual flightpath, has a steady-state gain 
K-j = .190. 

Figure VI-6 shows the optimal gain K-j(k) and the adaptive filter 
gains K-j(k) for the simulation. No plot is shown for Magi 1 1 1 s algorithm, 
which computes the adaptive state estimate from (V-19) as a. weighted sum 
of parallel filter estimates and consequently does not use K^(k). -The 
optimal gain goes through three near-step changes; at t = 55 sec., 115 
sec., and 200 sec. For all adaptive gain plots shown, K^(k) lags K^(k) 
at each of these three times in changing to the new gain level. This 
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time lag is most noticeable in the filter of Sage and Husa, which has the 
most difficulty in estimating the optimal gain. We recall that the adap- 
tive filters all use sample innovations covariances, either of an assumed 
optimal filter in the method of Sage and Husa, or of each of a* group of 
parallel filters for the other methods. The true innovations covariances 
change immediately whenever K-j(k) changes, but the sample covariances are 
time averages and consequently change more slowly. The time lag in K-| (k) 
is most critical at t'= 55 sec. Here the adaptive gain remains low when 

the optimal gain is high, a condition which can produce high mean square 
* 

A • 

errors in e(k[k) (see Figure 1 1 1 - 1 ) . 

Figure VI-7 shows the error e(k) - y(k) in the envelope processor, 
as well as the error 0 {k) - e(kjk) for the optimal and constant-Q fil- 
ters. The reduction in error produced by the optimal filter is obvious. 
The nonadaptive filter with gain = .476 (Q set to Q^^,) reduces the 
error in y(k), but not as well as the optimal filter. The filter with 
gain of J 90 works about as well as the optimal, except for t > 200 
seconds, where Q and K^(k) go to zero. Figure VI-8 compares the optimal 
filter's error in ©{ k j k ) with that of the adaptive filters. All of these 
filters seem to work about as well as the optimal, except near the end of 
the approach, when Q(k) goes to zero. - 

We realize that results of a single simulation run cannot provide 
us with a firm basis for any meaningful conclusions. The results shown 
here are to some degree dependent upon the noise sequences tu(k) and v ( k ) 
peculiar to this particular run. We therefore repeat the above simula- 
tion 100 times: the simulation is repeated without reinitializing the 

machine random number generator, so that noise samples used in one 
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experimental run are independent of those used in the other runs. We 

A 

then obtain an ensemble average of the mean square error in ©( k | k ) for a 
given candidate filter. Let 0 ^ (kjk) be the error in ©( k j k) for the ith 
simulation run: 

( k ] k } = e(k) - e(kjk); ith run (VI-11) 

We obtain a sample mean square error by averaging the square error at 
time k for all simulation runs: 

i 100- . z . 

P( k l k) = —■ e ej(kjk) (VI-12) 

i=l 

Assuming e(kjk) to be unbiased, P(kjk) is a variance estimate. From (V- 
4) we know that the standard deviation in this estimate is given by 

a p = P(k|k)//T0072' = . 1 41 P( k I k) (VI-13) 

where P ( k | k) is the true error covariance for ©( k { k ) (which we know only 
for the optimal filter). We thus expect our sample error covariance 
P ( k j k ) to be within 14 percent of the true covariance most of the time. 

A A 

Here we use the square root of P ( k | k ) as a sample r.m.s. error in Q(k|k). 

A 

Figure VI-9 shows the sample r.m.s. error in 0 (k|k) computed from 
100 simulation runs for each of the estimators tested. The error in the 
envelope processor estimate y(k) stays near . 01 °; this is consistent with 
a constant covariance R(k) of 10**^. The optimal filter significantly 

lowers the error, usually holding it below .004°. The minimum innova- 

/ 

tions covariance and Alspach adaptive filters work nearly as well as the 
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optimal; only for t > 200 seconds, where Q( k) goes to zero, does the 
optimal filter have a noticeably smaller r.m.s. error. The adaptive 
filter of Sage-Husa does almost as well, with a peak r.m.s. error of 
about .006°. The Magi 1 1 filter was not tested here. It was much slower 
than the other adaptive filters, and the time required to run a 100- 
record simulation was too great to be practical (recall that the bank of 
parallel stationary filters must be implemented serially in simulation). 
We note, however, that the single-run results for this filter look very 
much like those of Alspach. This is not surprising, as both methods com- 
pute the a posteriori gain density p[K^ | from. parallel filters. 

A 

We note in Figure VI-9 that the constant-Q filter with gain K-j of 
.476 has an r.m.s. error of about .006°. This is somewhat higher than 
the error associated with the adaptive filters. The nonadaptive filter 

A 

with K-j at .190 does very well', however, with the r.m.s. error staying 
near .004°. This compares favorably with the optimal and adaptive fil- 
ters, Only when Q(k) becomes very small, as for t < 50 seconds and 
t > 200 seconds, do the adaptive filters work significantly better than 
this filter. 

Simulation Testing: The Deterministic Case 

We have established that the candidate adaptive and constant-Q fil- 
ters work reasonably well when the assumed state model is implemented. 

We must now find out how well they can work in a true physical environ- 
ment, where the aircraft is actually moving along a given flightpath. 

This of course is our original objective; to find an adaptive Kalman 
filtering scheme for computing a minimum mean square error estimate 
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0 ( k ] k ) , using the estimates output by the airborne receiver's envelope 
processor. 

The entire simulation of the preceding section has been repeated, 
with the logical variable MODEL now changed to FALSE. g( k ) is now 
updated deterministically as the aircraft travels at 120 knots along the 

y 

S-curve approach of Figure VI-1 (Recall that ©(t), ©(t), and ©(t) are 

shown for this flightpath in Figure VI-4). R(k) is still held constant 

-4 

at 10 , giving an r.m.s. error of .01° in the envelope processor esti- 

mate y(k). In addition to the candidate adaptive and constant-Q filters 
we also run the same optimal filter as before, setting Q(k - 1) to 
At8(tjJ. Of course, this filter is unrealizable, since ©(t^) is unknown 
to the aircraft. 

Results for the single simulation run are given in Figures VI-10, 
11, T2. Figure VI- TO depicts the gain K,(k) for the optimal and adaptive 
Kalman filters. The adaptive gains look about the same as for the sto- 
chastic model simulation, with the Sage-Husa filter again having the most 
difficulty in estimating the optimal gain. 

Figure VI-11 gives the error in ©(kjk) for the envelope processor 
as well as the optimal and constant-Q filters. Again, these plots look 
about the same as for the stochastic case. The optimal filter signifi- 
cantly lowers the envelope processor error, while the filter with con- 
stant gain at-. 476 also lowers the original error, but not as much. 

The filter with K-j at .190 performs about as well as the optimal . filter, 
except for t > 200 seconds. Here, where the aircraft is on runway cen- 
terline, the optimal filter has less error. Figure VI-12 compares the 
optimal filter's error with that of the adaptive filters. The minimum 
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innovations covariance, Alspach, and Magi 11 filters compare favorably 
with the optimal. The Sage-Husa filter experiences some difficulty, how- 
ever, having a noticeable error bias, especially in the time region 
between 50 and 100 seconds. 

Figure VI-13 shows the sample r.m.s. error in 0(k|k) .for 100 simu- 
lation runs for each of the filters tested. The r.m.s. error in the 
envelope processor estimate y(k) stays near .01°, as expected. The mini- 
mum innovations covariance and Alspach filters significantly reduce this 
error, generally holding it to .004° or less. Both filters show a brief 
rise in error at t = 55 seconds, where the acceleration ©(t) suddenly 
changes' (see Figure VI-4.C). Here the adaptive gain K^(k) lags the opti- 
mal gain, using a low suboptimal gain until the sample innovations covar- 
iances in the parallel filter bank can respond to the change in 6{t). 
Except for this temporary error increase, these two filters work about as 
well as the assumed optimal filter. The adaptive algorithm of Sage-Husa 
lowers the r.m.s. error in the envelope processor, but not to the same 
degree as the other adaptive filters. The error increase at t = 55 sec- 
onds is much more pronounced, rising almost to .012°. 


The nonadaptive filter with K-j at .476 has an r.m.s. error of 
about .006°. This is higher than the .004° error often present with the 


best adaptive filters. A reduction of error from .006° to .004° seems 


rather marginal, though, when we consider the added sophistication 
required by the adaptive filters. If we assume that maximum ©(t) is 

A 

known for the actual flightpath, then the constant-Q filter with K-j at 
.190 can be used. This filter has an r.m.s. error of about .004°; this 


error performance is only slightly different from that of the adaptive 




A. Envelope Processor; B. Optimal; C. Minimum Innov. Cov.'; D. Alspach 
E. Sage - Husa; F. = .476; G.- = .190 


Figure VI - 13. Sample KMS Error in 6(klk); 6{k)set by S- curve flightpath. 
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filter's. We note that the adaptive filters lower the r.m.s. 'error to 
about .002° when t is less than 50 seconds or greater than. 200 seconds. 
Citing Figure VI-4.C, e(t) is nearly zero for t < 50 seconds, and it 
equals zero for t > 200 seconds (here the aircraft is flying ; . on runway 
centerline). Thus the adaptive filters work best only when ithe -accel era- 
tion is near zero. 

From these results it would seem wise to use the constant-Q filter 
with knowledge of maximum e(t) for estimating © ( k) . This filter works 
nearly as well as the best adaptive filters, and is much simpler to 
- implement. The constant-Q filter uses the same equations as the Kalman 
filter of- (II 1-13) -(II 1-28) . On the other hand, the minimum innovations 
covariance and Alspach adaptive filters, which have the best error per- 
formance, must impiement parallel Kalman filters. While such parallel 
processing would be fast, especially for the minimum innovations covari- 
ance filter, the hardware cost involved In realizing a bank of parallel 
filters seems unjustified by the marginal improvement in estimation 
error. 

• A 

The constant-Q filter use here with gain K-j of .190 requires a 
knowledge of the maximum acceleration for the actual landing 
approach. In the general problem statement: of this'paper we have assumed 
the landing approach to be unknown;, we can compute on ^ f° r 
family of allowable flightpaths. We therefore cannot realize this filter 
for our estimation problem as formally described. Yet we assume that the 
approach pattern used at a given airport and runway is standard. In such 
a case could be computed and stored until needed by an aircraft 

landing at that runway. At the beginning of the standard approach, this 



value e MAX could be transmitted to the aircraft, which could then, use the 
above constant-Q filter. . 

If the landing approach is not standard, or if there is no provi- 
sion for communicating © MA v to the aircraft, we can still use the 
constant-Q filter which uses 0^^ for the family of allowable flightpaths 
(the filter with gain K.| = .476 in our simulation). This filter has a 
higher r.m.s. error in e(k|k) than the best adaptive filters: ,006° com- 

pared to .004° when the aircraft is maneuvering. Yet the improvement in 
performance for, the adaptive filters is still not very significant when 
weighed against their added complexity. 

. In' general the adaptive filters work well. The filter of Alspach 
and the minimum innovations covariance filter lower the envelope proces- 
sor error in estimating s(k) from .01° to .004° or less: a 60% reduc- 

tion. But for our specific problem adaptive filtering does not appear to 
be necessary. The constant-Q filters work nearly as well and greatly 
simplify the estimation procedure. We should point out that use of the 
constant-Q filters is made possible by the fact that the measurement 
noise covariance R( k) is known in our problem. In the general adaptive 
estimation problem, where Q(k) and R(k) are both unknown, we would expect 
the constant-Q filter to be highly suboptimal , with the superiority of 
the adaptive Kalman filters becoming clearly evident. We note that the 
minimum innovations covariance and Alspach filters do not use R(k) 
anyway. These filters would work just as well for our problem if R(k) 
had been unknown. 

We digress here to make a general observation regarding adaptive 
filtering which may be useful in future work. Both in the simulation 
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runs- presented here and in other runs the filter of Sage and Husa was 

i . ^ 

slower' than the other adaptive filters in responding to changes in Q(k). 
The adaptive gain K(k) did not follow the optimal gain K(k) as well, and 

. V A 

the time lags between K{k) and K(k) were much more pronounced when Q(k) 
varied, .rapidly. We can speculate as to why this happened. .-The Sage-Husa 
filter uses the sample innovations covariance of the adaptive filter in 
computing K(k). If Q(k) experiences a step change, the sample covariance 

A 

eventually detects this change, causing K(k) to move toward the new 
steady-state gain. Yet the sample innovations covariance cannot immedi- 

A 

ately show the effect of using the new adaptive gain K(k), since most of 
the innovations residuals used in this statistic are those computed when 

At 

K(k) was at the old value. The other adaptive algorithms use banks of 
fixed-gain parallel Kalman filters. When the adaptive gain moves to a 

A 

new value K(k) as the result of a change -in Q(k), we already have a par- 
allel filter operating with a'gain close to K(k). This parallel filter 

. , A 

has always been running at the same gain., so that, the effects of K(k) 
upon the sample innovations covariance for the new value of Q ( k) are felt 
much sooner. 

We should bear in mind that the filter of Sage-Husa was designed 
for problems where the noise was stationary. Also,, we elected to use the 
suboptimal algorithm of Sage-Husa, rather than their more complicated 
optimal design based on MAP estimation. We speculate that this filter, 
by using smoothed state estimates, would be much faster in adapting to 
changes in Q(k) and R(k). 



CHAPTER VII 


CONCLUSION 

We have examined adaptive Kalman filtering for use in estimating an 
aircraft's azimuth angle e(k) in the Microwave Landing System. Adaptive 
filters from the literature were modified for application to the MLS 
problem and then tested in a simulated landing approach. The airborne 
receiver's envelope processor azimuth estimate y(k) was used as an input 
to each candidate filter. The filter's task was to produce a new esti- 

A 

mate e(k|k) having less mean-square error than y(k). In the simulation 
testing conducted here, where an S-curve flightpath was used, two adap- 
tive filters performed well: the r.m.s. error in y(k) was. lowered during 

various pnases of the approach by 60 percent or more. A suboptimal , non- 
adaptive estimation scheme was found to work almost as well. 

In Figure VI-13 we presented the results of our S-curve landing 
simulation, where the square error in e(kjk) was averaged for time k over 
100 simulated approaches. The minimum innovations covariance filter and 
the filter of Alspach proved to be the best adaptive filters, generally 
holding the r.m.s. error in ©( k j k ) to .004° or less. This compares to a 
constant r.m.s. error in the envelope processor estimate y(k) of ,01°. 

An r.m.s. error of about .004°' was obtained by a suboptimal filter using 
a fixed estimate of the state noise covariance Q(k) in the assumed sto- 
chastic model for 0 (k). The estimate of Q was based upon knowledge of 
the maximum acceleration © MAX for the actual S-curve flightpath. When 

M •* 

for the flightpath was unknown a second value of 0,, nv was used, 

MAX 


115 
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based on maximum acceleration for an allowable family of flightpaths. A 
filter using a constant $ estimate based on this value of 0 ^^ estimated - 
e(k) with an r.m.s. error of about .006°. 

In Figure VI-13 we note that the adaptive filters lower the r.m.s. 

A < 

error to about .002° for t above 230 seconds; here the aircraft has been 
on runway centerline for about a mile. This is the only phase of the 
landing approach where the adaptive filters significantly outperform the 
best fixed-Q filter (the filter with steady-state gain K-j of .190 for our 
simulation). 

While this constant-Q filter has an error performance comparable to 
that of the adaptive filters, it is not realizable under the formal con- 
straints of our 1 problem as defined in Chapter II. We assumed that the 
flightpath of the aircraft was unknown to the candidate filter; only the 
restrictions on the family of allowable flightpaths were given. Thus we 
would not know Sjv^v for the actual landing approach. As explained in 
Chapter VI, however, this filter could be used at an airport runway where 
the landing approach is standard. © MAX could be computed and stored for 
a given standard approach, and its value subsequently transmitted to an 
approaching aircraft. 

If there is no provision for making the value 0 ,^ tor the given 
flightpath available to the aircraft's MLS receiver, we can resort to the 


less optimal constant-Q filter where 0 ^^ for the set of allowable 
flightpaths is used.. In our simulation this filter compared reasonably 
well with the best adaptive filters in error performance (.006° r.m.s. 
error vs. .004° r.m.s. error for large segments of the approach). 

The fixed-Q filters use the Kalman filter equations of (III-18)- 



117 


(III -28) , replacing Q(k) with Q,^. Revisiting the algorithm flowcharts 
in Figures V-l to V-4, we note that the adaptive Kalman filters are much 
more complex. The minimum innovations covariance and Alspach filters can 
run nearly as fast as the Kalman filter, but only when the parallel sta- 
tionary filters run simultaneously. This requires the use of a special- 
purpose digital machine with parallel processing capability. 

We conclude that, for a curved landing pattern similar to the S- 
curve approach used here, adaptive Kalman filtering i.s not needed. A 
nonadaptive, fixed-Q filter can estimate 0(k) with a mean-square error 
performance comparable to optimal. While adaptive filters can lower this 
mean-square error further, the marginal improvement is judged to be 
insignificant in comparison with the added complexity and cost of realiz- 
ing a bank of parallel filters. 

The success of the' constant-Q filter results from the fact that the 
measurement noise covariance R(k) is known in our problem. Given a know- 
ledge of R{k), we can be assured that. our nonadaptive filter will not 
diverge by merely setting our estimate of Q(k) to some maximum limit 
never exceeded by the true value. The more accurately we can bound Q(k), 
the closer to optimal this filter becomes. In the more general problem 
where Q(k) .and R(k) are both unknown, we speculate that an adaptive 
filter would be clearly superior. We note that the minimum innovations 
covariance and Alspach adaptive filters make no use of R ( k } in our 
problem, but assume it to be unknown. Thus their error performance would 
remain unchanged if knowledge of R(k) were lost. 
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